
# CausalANOVAFit ! adjusted -------------------------------------------------------------

### ADJUSTED VERSION OF FUNCTION

CausalANOVAFit <- function(formula, internal.int,
                           data,cost,pair.id=NULL,
                           nway=2, diff=TRUE, eps=1e-5, family="binomial",
                           fac.level=NULL,ord.fac=NULL, cluster=cluster,
                           verbose=TRUE,
                           screen=FALSE, screen.type="fixed", screen.num.int=5,
                           collapse=FALSE, collapse.type="fixed", collapse.cost=0.3,
                           maxIter=100){
    
    ## ############################
    ## (0.1) Setup
    ## ############################ 
    Gorder <- nway
    formula.orig <- formula
    ## We convert the data into the data only with necessary variables. 
    data.main <- model.frame(formula, data=data)
    data.x <- model.frame(formula, data=data)[,-1]
    
    ## Extract information.
    if(is.null(fac.level)){
        fac.level <- unlist(lapply(data.x,FUN=function(x) length(levels(x))))
    }
    n.fac <- length(fac.level)
    
    if(is.null(ord.fac)){
        ord.fac <- unlist(lapply(data.x,FUN=function(x) is.ordered(x)))
    }
    for(i in 1:ncol(data.x)){
        if(ord.fac[i]==TRUE){
            data.main[,(i+1)] <- factor(data.main[,(i+1)],ordered=FALSE)
        }
    }
    
    if(any(ord.fac[fac.level==2]==TRUE)){
        if(verbose==TRUE) cat("Note: binary factors are recoded to unordered.")
        ord.fac[fac.level==2] <- FALSE
    }
    
    ## Print the status
    if(verbose==TRUE){
        print.fac <- cbind(fac.level,as.data.frame(ord.fac))
        colnames(print.fac) <- c("levels","ordered")
        cat("\nCheck: the number of levels for factors and whether they are ordered.\n")
        print(print.fac);rm(print.fac)
    }
    
    
    
    ## ############################
    ## (0.2) Screening
    ## ############################
    if(Gorder!=1){
        screen.out <- ScreenINT(formula=formula, data.main=data.main,
                                data.x=data.x, screen=screen, screen.type=screen.type,
                                family=family, fac.level=fac.level, screen.num.int=screen.num.int,
                                Gorder=Gorder, maxIter=maxIter, internal.int=internal.int,
                                verbose=verbose)
        
        indTwo <-  screen.out$indTwo
        indThree <- screen.out$indThree
        formula <-  screen.out$formula
        Gorder <- screen.out$Gorder
        
        if(screen==FALSE){
            formula.orig <- formula
        }
    }else{
        indTwo <-  NULL
        indThree <- NULL
    }
    
    rm(data.x); 
    data <- model.frame(formula,data=data.main); rm(data.main)
    level.list <- lapply(data[,-1], levels)
    data.orig <- data
    order.f <- attr(terms(formula,data=data), "order")
    
    ## ############################
    ## (0.3) Take Differences for Paired Data
    ## ############################ 
    if(diff==TRUE){
        data <- data[order(pair.id),]
        side <- rep(c(1,0),times=nrow(data)/2)
    }else{
        side <- NULL
    }
    
    if(diff==TRUE){
        data1 <- data[side==1,]
        data2 <- data[side==0,]
        cluster <- cluster[side==1]
        X1 <- model.matrixBayes(formula,data=data1)
        X2 <- model.matrixBayes(formula,data=data2)
        X <- X1 - X2 
        y <- model.frame(formula,data=data1)[,1]
        base.name <- colnames(X1)
        if(cost==1){
            model.frameC  <- model.frame(formula,data=data1)
            contr <- rep(list("contr.sum"), ncol(model.frameC) - 1)
            names(contr) <- colnames(model.frameC)[-1]
            x1 <- model.matrix(formula, data=data1,contrast=contr)[,-1]
            x2 <- model.matrix(formula, data=data2,contrast=contr)[,-1]
            X.no <- x1-x2; rm(model.frameC)
        }
    }else if(diff==FALSE){
        X <- model.matrixBayes(formula,data=data)
        y <- model.frame(formula,data=data)[,1]
        base.name <- colnames(X)
        if(cost==1){
            model.frameC  <- model.frame(formula,data=data)
            contr <- rep(list("contr.sum"), ncol(model.frameC) - 1)
            names(contr) <- colnames(model.frameC)[-1]
            X.no <- model.matrix(formula, data=data,contrast=contr)[,-1]
            rm(model.frameC)
        }
    }
    
    ## ############################ 
    ## (0.4) If No Regularization
    ## ############################ 
    if(cost==1 & screen==FALSE){
        ## if(verbose==TRUE) cat("\nNo Collapsing when 'collapse=FALSE'.\n")      
        NoReg <- NoRegularization(y=y, X.no=X.no, base.name=base.name, fac.level=fac.level,
                                  ord.fac=ord.fac, Gorder=Gorder, indTwo=indTwo,indThree=indThree,
                                  cluster=cluster)                
        intercept <- NoReg$intercept
        coefs <- NoReg$coefs
        vcov  <- NoReg$vcov
        
        ## Index
        levelIndex <- CreatelevelIndex(fac.level=fac.level, ord.fac=ord.fac, Gorder=Gorder,
                                       indTwo=indTwo, indThree=indThree)
        use.ind <-  (levelIndex$plus==1)*(levelIndex$dif==0)
        Index.use <- levelIndex[use.ind==1,]
        Index.use$start <- c(1,cumsum(Index.use$length)[-nrow(Index.use)]+1)
        Index.use$end <-  cumsum(Index.use$length)
        
        order.f <- attr(terms(formula, data=data), "order")
        
        ###### ATTEMPT AT ADAPTED VERSION
        CI.tab <- list()
        if(is.null(refcat)){
            refcat <- rep(1, length(coefs))
        }
        if(length(refcat) != length(coefs)){
            stop("Length of reference categories incorrect")
        }
        for(z in 1:length(coefs)){
            seq.u <- Index.use$start[z]:Index.use$end[z]
            effect <- coefs[[z]] - coefs[[z]][refcat[z]]
            vcov.p <- c()
            for(i in 1:(length(seq.u))){
                vcov.p[i] <- VarEffect(ind1=seq.u[i],ind2=seq.u[refcat[z]],vcov.full=vcov)
            }
            sd <- sqrt(round(vcov.p,digits=10))
            table <- as.data.frame(cbind(effect, sd, effect - 1.96*sd, effect + 1.96*sd))
            if(Index.use$order[z]==1){
                colnames(table) <- c("AME", "sd", "2.5% CI", "97.5% CI")
            }else{
                colnames(table) <- c("AMIE", "sd", "2.5% CI", "97.5% CI")
            }
            CI.tab[[z]] <- table
        }
        
        AME <- coefs[order.f==1]
        AME <- lapply(AME, function(x) round(x/(eps*10))*(eps*10))
        if(Gorder>=2){
            AMIE2 <- coefs[order.f==2]
            AMIE2 <- lapply(AMIE2, function(x) round(x/(eps*10))*(eps*10))
        }else{
            AMIE2 <- NULL
        }
        if(Gorder==3){
            AMIE3 <- coefs[order.f==3]
            AMIE3 <- lapply(AMIE3, function(x) round(x/(eps*10))*(eps*10))
        }else{
            AMIE3 <- NULL
        }
        
        output <- list("intercept"=intercept,
                       "coefs"=coefs,"data"=data,
                       "vcov"=vcov, "CI.table"=CI.tab,
                       "AME"=AME, "AMIE2"=AMIE2,"AMIE3"=AMIE3,
                       "nway"=nway,"formula"=formula,
                       "cost"=cost,"fac.level"=fac.level,
                       "ord.fac"=ord.fac,
                       "level.list"=level.list,
                       "diff"=diff,
                       "eps"=eps,"side"=side,
                       "data.orig"=data.orig,
                       "formula.orig"=formula.orig,
                       "pair.id"=pair.id,
                       "indTwo"=indTwo,
                       "indThree"=indThree,
                       "family"=family,
                       "internal.int"=internal.int,
                       "Gorder"=Gorder,
                       "inference"=TRUE,
                       "screen"=screen,
                       "screen.type"=screen.type,
                       "screen.num.int"=screen.num.int,
                       "collapse"=collapse,
                       "collapse.type"=collapse.type,
                       "collapse.cost"=collapse.cost)
        
        class(output) <- c("CausalANOVAFit","list")
        return(output)
    }
    
    ## ##################
    ## (1.0) Create the basic Index
    ## ################## 
    levelIndex <- CreatelevelIndex(fac.level=fac.level,ord.fac=ord.fac, 
                                   Gorder=Gorder, indTwo=indTwo, indThree=indThree)   
    coef.length <- sum(levelIndex$length)
    coef.length.mu <- coef.length + 1
    
    
    ## #########################
    ## Making the basic matrices
    ## Note: these do not change depending on X. Only depends on fac.level, ord.fac, indTwo, indThree
    ## ######################### 
    ## (1.1) Making L matrix    
    L <- Lcombinefunction(fac.level=fac.level,ord.fac=ord.fac,facCons=FALSE,Gorder=Gorder,
                          indTwo=indTwo,indThree=indThree)
    
    
    ## (1.2) Making Psy matrix
    PsyConstraintMat <- PsyConstraintCombine(fac.level=fac.level,ord.fac=ord.fac,
                                             Gorder=Gorder, 
                                             indTwo=indTwo, indThree=indThree)
    
    ## (1.3): Making the weights
    weight.fac <- c()
    weight.ols <- c()
    for(w in 1:n.fac){
        weight.list <- CreateWeights(formula,data=data,Gorder=Gorder,
                                     facCons=FALSE,
                                     side=side,fac.level=fac.level,dif=diff,
                                     ord.fac=ord.fac,type.ind=w,
                                     indTwo=indTwo, indThree=indThree)
        weight.fac <- c(weight.fac, weight.list$weight.fac)
        weight.ols <- c(weight.ols, weight.list$weight.ols)        
    }
    ## weight.fac <- (weight.fac/sum(weight.fac))*length(weight.fac)
    ## weight.fac <- (weight.fac/sum(weight.fac))*length(weight.fac)
    weight <- weight.fac * weight.ols
    
    ## (1.4): Create Slack variables
    l.slack <- lengthSlack(fac.level=fac.level,ord.fac=ord.fac,Gorder=Gorder,
                           facCons=FALSE)$l.slack
    
    ## (1.5) : Normalize sum constraints so that 0<t<1
    sum.constraint <- c(rep(0,times=coef.length.mu), weight)
    t.sum <- cost*sum(weight.fac)
    
    ## (1.6): Positivity constraint
    Positive <- diag(coef.length.mu+l.slack)
    
    ## (1.7): Create zero-sum constraint
    Constraint <- CreateANOVAconst(fac.level=fac.level,ord.fac=ord.fac,
                                   Gorder=Gorder,facCons=FALSE,
                                   indTwo=indTwo, indThree=indThree)
    
    ## (1.8) : Create Z matrix (Expand version of Design Matrix)
    Zcombine <- Zcombinefunction(X=X,fac.level=fac.level,ord.fac=ord.fac, 
                                 Gorder=Gorder,
                                 indTwo=indTwo, indThree=indThree)
    ## Adjust Z matrix for slack variables
    ZcombineF <- cbind(Zcombine, matrix(0,ncol=l.slack,nrow=nrow(Zcombine)))
    
    ## (1.9) : Combine all inequality constraints    
    G <- rbind(-PsyConstraintMat, -sum.constraint, Positive)
    H <- c(rep(0,times=nrow(PsyConstraintMat)),-t.sum,rep(0,times=nrow(Positive)))
    LargeG <- rbind(L, -L, Constraint, -Constraint, G)
    LargeH <- c(rep(-eps, times=nrow(L)),rep(-eps, times=nrow(L)),
                rep(-eps, times=nrow(Constraint)),rep(-eps, times=nrow(Constraint)),
                H)
    
    if(collapse==TRUE){if(verbose==TRUE) cat("\n ***** Collapsing ***** \n")}
    
    ## ############################
    ## (2.0) Fit the Main Model
    ## ############################ 
    solve.fit <- Glsei(A=ZcombineF,B=y,G=LargeG,H=LargeH,verbose=TRUE,
                       tol=5*eps,type="2")
    
    BaseCoef <- solve.fit$X[-1]
    coefs <- CoefExtract(BaseCoef,base.name=base.name, 
                         fac.level=fac.level, ord.fac=ord.fac,
                         Gorder=Gorder,indTwo=indTwo,indThree=indThree)
    
    AME <- coefs[order.f==1]
    AME <- lapply(AME, function(x) round(x/(eps*10))*(eps*10))
    if(Gorder>=2){
        AMIE2 <- coefs[order.f==2]
        AMIE2 <- lapply(AMIE2, function(x) round(x/(eps*10))*(eps*10))
    }else{
        AMIE2 <- NULL
    }
    if(Gorder==3){
        AMIE3 <- coefs[order.f==3]
        AMIE3 <- lapply(AMIE3, function(x) round(x/(eps*10))*(eps*10))
    }else{
        AMIE3 <- NULL
    }
    
    output <- list("solve.fit"=solve.fit,
                   "intercept"=solve.fit$X[1],
                   "coefs"=coefs,"data"=data,
                   "AME"=AME, "AMIE2"=AMIE2,"AMIE3"=AMIE3,
                   "nway"=nway,"formula"=formula,
                   "cost"=cost,"fac.level"=fac.level,
                   "level.list"=level.list,
                   "ord.fac"=ord.fac,"diff"=diff,
                   "eps"=eps,"side"=side,
                   "data.orig"=data.orig,
                   "formula.orig"=formula.orig,
                   "pair.id"=pair.id,
                   "indTwo"=indTwo,
                   "indThree"=indThree,
                   "family"=family,
                   "screen"=screen,
                   "internal.int"=internal.int,
                   "screen.type"=screen.type,
                   "screen.num.int"=screen.num.int,
                   "Gorder"=Gorder,
                   "inference"=FALSE,
                   "collapse"=collapse,
                   "collapse.type"=collapse.type,
                   "collapse.cost"=collapse.cost)
    
    class(output) <- c("CausalANOVAFit","list")
    return(output)
}


# plotCausalANOVA !adjusted ---------------------------------------------------------

## ADJUSTED VERSION OF FUNCTION

plot.CausalANOVA <- function(x, fac.name,
                             treat.ind=1,
                             type = "ConditionalEffect", space = 15, xlim, 
                             maintitle = "AMIE", ylabs = rev(Names), col = "black",
                             ...){
    treat.ind <-  1
    adj.p   <- 2.3
    add     <- FALSE
    ## HouseKeeping 
    if(x$Gorder==1 | missing(fac.name)==TRUE){
        type <- "AME"
    }
    
    if(type=="AME"){
        formula <- x$formula
        terms.f <- terms(formula, data=x$data)
        order.f <- attr(terms.f, "order")
        var.name <- attr(terms.f, "term.labels")[order.f==1]
        level.list <- x$level.list        
        inference <- x$inference
        
        if(inference==TRUE){
            AME.CI  <- x$CI.table[which(order.f==1)]
            font <- CI.value <- CI.up <- CI.low <- c()
            Name <- list()
            for(z in 1:length(AME.CI)){
                Name0 <- c(paste("Factor: ", var.name[z], sep=""),
                           paste("    ", level.list[[z]],sep=""))
                Name <- c(Name, Name0)
                font <- c(font, c(2, rep(1, length(level.list[[z]]))))
                AME.base <- AME.CI[[z]]
                AME.base[,3] <- AME.base[,1] -1.96*AME.base[,2]
                AME.base[,4] <- AME.base[,1] +1.96*AME.base[,2]
                CI.value <- c(CI.value, c(NA, AME.base[,1]))
                CI.up <- c(CI.up, c(NA, AME.base[,4]))
                CI.low <- c(CI.low, c(NA, AME.base[,3]))               
            }
            xlim.min <- min(na.omit(CI.low))
            xlim.max <- max(na.omit(CI.up))
        }else{
            AME.CI  <- x$AME
            font <- CI.value <- c()
            Name <- list()
            for(z in 1:length(x$AME)){
                Name0 <- c(paste("Factor: ", var.name[z], sep=""),
                           paste("    ", level.list[[z]],sep=""))
                Name <- c(Name, Name0)
                font <- c(font, c(2, rep(1, length(level.list[[z]]))))
                ame <- unlist(AME.CI[[z]])
                ame <- ame - ame[length(ame)]
                CI.value <- c(CI.value, NA, ame)
            }
            xlim.min <- min(na.omit(CI.value))
            xlim.max <- max(na.omit(CI.value))
        }
        
        if(missing(xlim)==FALSE){
            par(mar=c(4,space,2,2))
            plot(rev(CI.value), seq(from=1, to=length(CI.value)), pch=19, yaxt="n", xlim=xlim,
                 ylim=c(0, length(CI.value)+1), ylab="", main="AME", xlab="Estimated Effects", col=col)
            if(inference==TRUE){
                segments(rev(CI.low), seq(from=1, to=length(CI.value)),
                         rev(CI.up), seq(from=1, to=length(CI.value)), lwd=2,col=col)
            }
            Axis(side=2, at=seq(from=1, to=length(CI.value)), labels=rep("",length(Name)))
            Axis(side=2, at=seq(from=1, to=length(CI.value)), labels=rev(Name),
                 las=2, hadj=0, line=10, tck=0, lwd=0)
            abline(v=0, lty=2)
        }else{
            par(mar=c(4,space,2,2))
            plot(rev(CI.value), seq(from=1, to=length(CI.value)), pch=19, yaxt="n",
                 xlim=c(xlim.min, xlim.max),
                 ylim=c(0, length(CI.value)+1), ylab="", main="AME", xlab="Estimated Effects", col=col)
            if(inference==TRUE){
                segments(rev(CI.low), seq(from=1, to=length(CI.value)),
                         rev(CI.up), seq(from=1, to=length(CI.value)), lwd=2,col=col)
            }
            Axis(side=2, at=seq(from=1, to=length(CI.value)), labels=rep("",length(Name)))
            Axis(side=2, at=seq(from=1, to=length(CI.value)), labels=rev(Name),
                 las=2, hadj=0, line=10, tck=0, lwd=0)
            abline(v=0, lty=2)
        }        
        
    }else if(type=="AMIE"){
        formula <- x$formula
        indTwo <- x$indTwo
        indThree <- x$indThree
        plot.Gorder <- length(fac.name)
        terms.f <- terms(formula,data=x$data)
        order.f <- attr(terms.f, "order")
        var.name <- attr(terms.f, "term.labels")[order.f==1]
        inference <- x$inference
        AMIE2 <- x$AMIE2
        AMIE3 <- x$AMIE3
        
        ## This accomodate Three-ways 
        Fac.ind <- c()
        for(z in 1:length(fac.name)){
            Fac.ind[z] <- which(fac.name[z]==var.name)
        }
        Fac.ind <- Fac.ind[order(Fac.ind)]
        
        ## Find the index
        if(plot.Gorder==2){
            INT.ind <- which(apply(c(Fac.ind[1], Fac.ind[2]) == indTwo, 2, all)==TRUE)
            if(inference==TRUE) AMIE.CI  <- x$CI.table[which(order.f==2)]
            else AMIE.CI  <- AMIE2
        }
        if(plot.Gorder==3){
            INT.ind <- which(apply(c(Fac.ind[1], Fac.ind[2], Fac.ind[3]) == indThree, 2, all)==TRUE)
            if(inference==TRUE) AMIE.CI  <- x$CI.table[which(order.f==3)]
            else AMIE.CI  <- AMIE3
        }
        
        ## Plot: AMIE (just need fac.name)
        if(inference==TRUE){
            CI.value <- AMIE.CI[[INT.ind]][,1]
            CI.up <- AMIE.CI[[INT.ind]][,4]
            CI.low <-  AMIE.CI[[INT.ind]][,3]
            Name   <- rownames(AMIE.CI[[INT.ind]])
            xlim.min <- min(na.omit(CI.low))
            xlim.max <- max(na.omit(CI.up))
        }else{
            amie <- unlist(AMIE.CI[[INT.ind]])
            CI.value <- amie - amie[length(amie)]
            Name <- names(CI.value)
            xlim.min <- min(na.omit(CI.value))
            xlim.max <- max(na.omit(CI.value))
        }
        
        
        # multiplying by 100 for percentage points
        CI.value <- CI.value * 100
        CI.low <- CI.low * 100
        CI.up <- CI.up * 100
        xlim.min <- xlim.min * 100
        xlim.max <- xlim.max * 100
        
        ## plot.value <- plot.value - plot.value[base.ind]
        
        max.nchar <- max(nchar(ylabs))/adj.p
        ## Create blank plot first
        if(missing(xlim)==FALSE){
            par(mar=c(4,space,2,2))
            plot(rev(CI.value), seq(from=1, to=length(CI.value)), pch=19, yaxt="n", xlim=xlim,
                 ylim=c(0, length(CI.value)+1), ylab="", main=maintitle, 
                 xlab="Change in Pr(Policy Chosen) in percentage points",col=col)
            if(inference==TRUE){
                segments(rev(CI.low), seq(from=1, to=length(CI.value)),
                         rev(CI.up), seq(from=1, to=length(CI.value)), lwd=2, col=col)
            }
            Axis(side=2, at=seq(from=1, to=length(CI.value)), labels=rep("",length(Name)))
            Axis(side=2, at=seq(from=1, to=length(CI.value)), labels=ylabs,
                 las=2, hadj=0, line=max.nchar, tck=0, lwd=0)
            abline(v=0, lty=2)                         
        }else{
            par(mar=c(4,space,2,2))
            plot(rev(CI.value), seq(from=1, to=length(CI.value)), pch=19, yaxt="n",
                 xlim=c(xlim.min, xlim.max),
                 ylim=c(0, length(CI.value)+1), ylab="", main=maintitle, 
                 xlab="Change in Pr(Policy Chosen) in percentage points",col=col)
            if(inference==TRUE){
                segments(rev(CI.low), seq(from=1, to=length(CI.value)),
                         rev(CI.up), seq(from=1, to=length(CI.value)), lwd=2, col=col)
            }
            Axis(side=2, at=seq(from=1, to=length(CI.value)), labels=rep("",length(Name)))
            Axis(side=2, at=seq(from=1, to=length(CI.value)), labels=ylabs,
                 las=2, hadj=0, line=max.nchar, tck=0, lwd=0)
            abline(v=0, lty=2)
        }
        
    }else if(type=="ConditionalEffect"){      
        treat.fac <- fac.name[treat.ind]
        if(treat.ind==2) cond.ind <- 1
        else if(treat.ind==1) cond.ind <- 2
        cond.fac  <- fac.name[cond.ind]
        inference <- x$inference
        CE <- ConditionalEffect(object=x, treat.fac=treat.fac, inference=inference,
                                cond.fac=cond.fac, base.ind=NULL, verbose=FALSE)
        CE.plot <- CE$Conditional
        cond.level <- CE$cond.level
        treat.level <- CE$treat.level
        
        if(inference==TRUE){
            font <- CE.value <- CE.up <- CE.low <- c()
            Name <- list()
            for(z in 1:length(CE.plot)){
                Name0 <- c(paste(cond.fac, "=", cond.level[z]),
                           paste("    ", treat.level,sep=""))
                Name <- c(Name, Name0)
                font <- c(font, c(2, rep(1, length(treat.level))))
                CE.value <- c(CE.value, c(NA, CE.plot[[z]][,1]))
                CE.up <- c(CE.up, c(NA, CE.plot[[z]][,4]))
                CE.low <- c(CE.low, c(NA, CE.plot[[z]][,3]))
            }
            xlim.min <- min(na.omit(CE.low))
            xlim.max <- max(na.omit(CE.up))
            
            if(missing(xlim)==FALSE){
                par(mar=c(4,space,2,2))
                plot(rev(CE.value), seq(from=1, to=length(CE.value)), pch=19, yaxt="n", xlim=xlim,
                     ylim=c(0, length(CE.value)+1), ylab="", xlab = "Estimated Effects", 
                     col=col, main = "Conditional Effects")
                segments(rev(CE.low), seq(from=1, to=length(CE.value)),
                         rev(CE.up), seq(from=1, to=length(CE.value)), lwd=2, col=col)
                Axis(side=2, at=seq(from=1, to=length(CE.value)), labels=rep("",length(Name)))
                Axis(side=2, at=seq(from=1, to=length(CE.value)), labels=rev(Name),
                     las=2, hadj=0, line=10, tck=0, lwd=0)
                abline(v=0, lty=2)
            }else{
                par(mar=c(4,space,2,2))
                plot(rev(CE.value), seq(from=1, to=length(CE.value)), pch=19, yaxt="n",
                     xlim=c(xlim.min, xlim.max),
                     ylim=c(0, length(CE.value)+1), ylab="", 
                     xlab="Estimated Effects", main = "Conditional Effects", col=col)
                segments(rev(CE.low), seq(from=1, to=length(CE.value)),
                         rev(CE.up), seq(from=1, to=length(CE.value)), lwd=2, col=col)
                Axis(side=2, at=seq(from=1, to=length(CE.value)), labels=rep("",length(Name)))
                Axis(side=2, at=seq(from=1, to=length(CE.value)), labels=rev(Name),
                     las=2, hadj=0, line=10, tck=0, lwd=0)
                abline(v=0, lty=2) 
            }
        }else{
            font <- CE.value <- c()
            Name <- list()
            for(z in 1:length(CE.plot)){
                Name0 <- c(paste(cond.fac, "=", cond.level[z]),
                           paste("    ", treat.level,sep=""))
                Name <- c(Name, Name0)
                font <- c(font, c(2, rep(1, length(treat.level))))
                CE.value <- c(CE.value, c(NA, CE.plot[[z]]))
            }
            xlim.min <- min(na.omit(CE.value))
            xlim.max <- max(na.omit(CE.value))
            
            if(missing(xlim)==FALSE){
                if(add==FALSE){
                    par(mar=c(4,space,2,2))
                }                    
                plot(rev(CE.value), seq(from=1, to=length(CE.value)), pch=19, yaxt="n", xlim=xlim,
                     ylim=c(0, length(CE.value)+1), ylab="", xlab="Estimated Effects", 
                     col=col, main = "Conditional Effects")
                Axis(side=2, at=seq(from=1, to=length(CE.value)), labels=rep("",length(Name)))
                Axis(side=2, at=seq(from=1, to=length(CE.value)), labels=rev(Name),
                     las=2, hadj=0, line=10, tck=0, lwd=0)
                abline(v=0, lty=2)
            }else{
                if(add==FALSE){
                    par(mar=c(4, space, 2,2))
                } 
                plot(rev(CE.value), seq(from=1, to=length(CE.value)), pch=19, yaxt="n",
                     xlim=c(xlim.min, xlim.max),
                     ylim=c(0, length(CE.value)+1), ylab="", 
                     xlab="Estimated Effects", col=col,
                     main = "Conditional Effects")
                Axis(side=2, at=seq(from=1, to=length(CE.value)), labels=rep("",length(Name)))
                Axis(side=2, at=seq(from=1, to=length(CE.value)), labels=rev(Name),
                     las=2, hadj=0, line=10, tck=0, lwd=0)
                abline(v=0, lty=2) 
            }
        }
    }    
}



# CausalANOVA -------------------------------------------------------------

CausalANOVA <- function(formula, int2.formula=NULL, int3.formula=NULL,
                        data, nway=1,
                        pair.id=NULL, diff=FALSE,
                        screen=FALSE, screen.type="fixed", screen.num.int=3,
                        collapse=FALSE, collapse.type="fixed", collapse.cost=0.3,
                        family="binomial", cluster=NULL,
                        maxIter=50, eps=1e-5,
                        fac.level=NULL, ord.fac=NULL,
                        select.prob=FALSE, boot=100, seed=1234,
                        verbose=TRUE){    
    
    ## House Keeping
    if(collapse==FALSE){
        cost <- 1
    }else if(collapse.type=="fixed"){
        cost <- collapse.cost
        if(cost > 1 | cost < 0 ){
            stop("Specify 'collapse.cost' between 0 and 1")
        }
    }
    if(missing(fac.level)) fac.level <- NULL
    if(missing(ord.fac)) ord.fac <- NULL        
    if(diff==TRUE & is.null(pair.id)==TRUE){
        stop("When 'diff=TRUE', specify 'pair.id'.")
    }
    if((nway %in% c(1, 2,3))==FALSE){
        stop("'nway' should be 1, 2 or 3.")
    }
    if(screen==TRUE & nway==1){
        warning("When 'nway=1', no screening is needed ('screen=FALSE').")
        screen <- FALSE
    }
    if((family %in% c("gaussian","binomial"))==FALSE){
        stop("'family' should be 'gaussian' or 'binomial'.")
    }
    y <- model.frame(formula,data=data)[,1]
    if(family=="gaussian" & is.numeric(y)==FALSE){
        stop("When 'family=gaussian', outcomes should be 'numeric'.") 
    }
    if(family=="binomial" & all(y %in% c(0,1))==FALSE){
        stop("When 'family=binomial', outcomes should be 0 or 1.") 
    }
    rm(y)
    if((screen.type %in% c("fixed","cv.min","cv.1Std"))==FALSE){
        stop("'screen.type' should be 'fixed', 'cv.min' or 'cv.1Std'.")
    }    
    if(is.null(int3.formula)==FALSE & nway!=3){
        warning("When 'int3.formula' is specified, 'nway=3'.")
        nway <- 3
    }
    if(is.null(int2.formula)==TRUE & is.null(int3.formula)==FALSE){
        stop("Cannot specify 'int3.formula' without specifying 'int2.formula'.")
    }
    if(is.null(int2.formula)==FALSE & screen==TRUE){
        warning("When 'int2.formula' is specified, 'screen' should be FALSE.")
        screen <- FALSE
    }
    if(nway==3 & screen==TRUE & screen.type=="fixed" & screen.num.int < 3){
        stop("'screen.num.int >=3' when 'nway=3'.")
    }    
    if(collapse.type=="cv.min" | collapse.type=="cv.1Std"){
        cat("\nRunning Cross Validation might take time...\n")
        cv <- cv.CausalANOVA(formula = formula,
                             int2.formula=int2.formula, int3.formula=int3.formula,
                             data=data,
                             cv.collapse.cost=c(0.1,0.3,0.7),
                             nfolds=5,
                             pair.id=pair.id, diff=diff,
                             nway=nway, family=family,
                             seed=seed, cluster=cluster,
                             screen=screen, screen.type=screen.type,
                             screen.num.int=screen.num.int,
                             maxIter=maxIter, eps=eps,                                                         
                             fac.level=fac.level,ord.fac=ord.fac, verbose=verbose)                            
        if(collapse.type=="cv.min") cost <- cv$cv.min
        if(collapse.type=="cv.1Std") cost <- cv$cv.1Std
        print(paste("Selected Cost parameter=", cost, sep=""))
    }
    if(cost==1 & select.prob==TRUE){
        select.prob <- FALSE
    }
    
    ## Factors only
    all.fac <- all(unlist(lapply(model.frame(formula,data=data)[,-1],
                                 FUN=function(x) is.element("factor",class(x)))))
    if(all.fac==FALSE) stop("Design matrix should contain only factors.")
    rm(all.fac);
    
    main.formula <- formula
    ## when int.formual != NULL
    if(is.null(int2.formula)==TRUE){
        internal.int <- TRUE
    }else if(is.null(int2.formula)==FALSE){
        internal.int <- FALSE
        if(is.null(int3.formula)==TRUE){
            formula <- as.formula(paste(as.character(formula)[2], "~",
                                        as.character(formula)[3], "+", as.character(int2.formula)[2], sep=""))
        }else if(is.null(int3.formula)==FALSE){
            formula <- as.formula(paste(as.character(formula)[2], "~",
                                        as.character(formula)[3],
                                        "+", as.character(int2.formula)[2],
                                        "+", as.character(int3.formula)[2],
                                        sep=""))
        }
    }
    
    ## ############################
    ## Main Function
    ## ############################ 
    fit <- CausalANOVAFit(formula=formula,
                          internal.int=internal.int,
                          data=data,cost=cost,pair.id=pair.id,
                          family=family, cluster=cluster,
                          screen=screen,screen.type=screen.type,
                          screen.num.int=screen.num.int,
                          maxIter=maxIter,
                          nway=nway,diff=diff,eps=eps,
                          collapse=collapse, collapse.type=collapse.type,
                          collapse.cost=collapse.cost,
                          fac.level=fac.level,ord.fac=ord.fac,
                          verbose=verbose)
    
    fit <- c(fit, main.formula=main.formula, int2.formula=int2.formula, int3.formula=int3.formula)
    
    if(select.prob==TRUE){
        stab.fit <- stab.CausalANOVA(fit,cluster=cluster,boot=boot,seed=seed)
        output <- list("fit"=fit, "stab.fit"=stab.fit)
        class(output) <- c("CausalANOVA","stab","list")
    }else{
        output <- fit
        class(output) <- c("CausalANOVA","list")
    }        
    return(output)
}


# CoefExtract -------------------------------------------------------------


CoefExtract <- function(BaseCoef, base.name, fac.level, ord.fac, Gorder, indTwo=NULL, indThree=NULL){
    ## input: BaseCoef
    ## I need to remove intercept first 
    
    n.fac <- length(fac.level)   
    levelIndex <- CreatelevelIndex(fac.level=fac.level,ord.fac=ord.fac,Gorder=Gorder,
                                   indTwo=indTwo, indThree=indThree)
    levelIndex$end <- cumsum(levelIndex$length)
    use.ind <- levelIndex$dif==0
    coefIndex <- levelIndex[use.ind==1,]
    plus.ind <- coefIndex$plus==1
    coefIndexPlus <- coefIndex[plus.ind==1,]
    coefIndexMinus <- coefIndex[plus.ind==0,]
    nameIndex <- coefIndexPlus
    nameIndex$start <- c(1, cumsum(nameIndex$length)[-nrow(nameIndex)]+1)
    nameIndex$end   <- cumsum(nameIndex$length)
    
    coefPlus <- c()
    Plus.list <- list()
    for(i in 1:nrow(coefIndexPlus)){
        ind <- coefIndexPlus$start[i]:coefIndexPlus$end[i]        
        Plus.list[[i]] <- BaseCoef[ind]       
        coefPlus <- c(coefPlus,BaseCoef[ind])
    }
    coefMinus <- c()
    Minus.list <- list()
    for(i in 1:nrow(coefIndexMinus)){
        ind <- coefIndexMinus$start[i]:coefIndexMinus$end[i]
        Minus.list[[i]] <- BaseCoef[ind]
        coefMinus <- c(coefMinus,BaseCoef[ind])
    }
    
    Coef.list <- list()
    for(i in 1:nrow(coefIndexPlus)){
        name.ind <- nameIndex$start[i]:nameIndex$end[i]
        coef <- Plus.list[[i]] - Minus.list[[i]]
        names(coef) <- base.name[name.ind]
        Coef.list[[i]] <- coef
    }    
    return(Coef.list)
}



# ConditionalEffect -------------------------------------------------------

ConditionalEffect <- function(object,treat.fac=NULL, cond.fac=NULL,
                              base.ind=1, round=3,
                              inference=NULL, verbose=TRUE){
    
    
    if(is.null(inference)==TRUE){
        inference <- object$inference
    }
    ## House Keeping
    ## only for Two-ways
    ## check mod
    ## check base.ind
    ## If selected interactions are not in indTwo.
    
    formula <- object$formula
    AME <- object$AME
    AMIE2 <- object$AMIE2
    fac.level <- object$fac.level
    ord.fac <- object$ord.fac
    Gorder <- object$Gorder
    indTwo <- object$indTwo
    indThree <- object$indThree
    data.main <- object$data[,-1]
    level.name <- lapply(data.main, levels)
    terms.f <- terms(formula,data=object$data)
    order.f <- attr(terms.f, "order")
    var.name <- attr(terms.f, "term.labels")[order.f==1]
    
    coefs.u <- unlist(object$coefs)
    
    if(inference==TRUE){
        vcov.u <- unlist(object$vcov)
    }
    
    ## LevelIndex
    ## Index
    levelIndex <- CreatelevelIndex(fac.level=fac.level, ord.fac=ord.fac, Gorder=Gorder,
                                   indTwo=indTwo, indThree=indThree)
    use.ind <-  (levelIndex$plus==1)*(levelIndex$dif==0)
    Index.use <- levelIndex[use.ind==1,]
    Index.use$start <- c(1,cumsum(Index.use$length)[-nrow(Index.use)]+1)
    Index.use$end <-  cumsum(Index.use$length)
    
    
    ## Find Index for Factors 
    treat.ind <- which(treat.fac==var.name)
    cond.ind <- which(cond.fac==var.name)
    treat.level <- levels(data.main[,treat.ind])
    cond.level <- levels(data.main[, cond.ind])
    rm(data.main)
    
    if(is.null(base.ind)==TRUE){
        base.ind <- length(treat.level)
    }
    
    Fac.ind <- c(cond.ind, treat.ind)
    Norotate <- all(order(Fac.ind) == c(1,2))
    Fac.ind <- Fac.ind[order(Fac.ind)]
    var.ind.mat <- matrix(seq(1:(fac.level[Fac.ind[1]]*fac.level[Fac.ind[2]])),
                          nrow=fac.level[Fac.ind[1]], ncol=fac.level[Fac.ind[2]])
    if(Norotate==FALSE){ var.ind.mat <- t(var.ind.mat) }
    
    
    ## Find the index for Main Effect
    AME.var.ind <- Index.use$start[treat.ind]:Index.use$end[treat.ind]
    AME.var.ind.mat <- as.matrix(cbind(AME.var.ind, AME.var.ind[base.ind]))
    ## rownames(AME.var.ind.mat) <- level.name[[treat.ind]]
    AME.var.ind.mat.f <- rep(list(AME.var.ind.mat), fac.level[cond.ind])
    ## names(AME.var.ind.mat.f) <- paste(var.name[cond.ind], "=", level.name[[cond.ind]],sep="")
    
    ## Find the index for Interaction 
    INT.ind <- which(apply(c(Fac.ind[1], Fac.ind[2]) == indTwo, 2, all)==TRUE)
    if(length(INT.ind)==0) stop("Specified Interactions are not in the model.")
    INT.ind.u <- sum(order.f==1) + INT.ind ## z
    AMIE.var.ind <- Index.use$start[INT.ind.u]:Index.use$end[INT.ind.u]
    
    ind.var <- list()
    for(i in 1:nrow(var.ind.mat)){
        ind.var0 <- matrix(NA, ncol=2, nrow=ncol(var.ind.mat))
        for(j in 1:ncol(var.ind.mat)){
            ind.var0[j, 1:2] <- c(var.ind.mat[i,j], var.ind.mat[i,base.ind])
        }
        ind.var[[i]] <- ind.var0
    }
    ## ## always row=Conditional, col=Treatment
    AMIE.var.ind.mat.f <- list()
    for(i in 1:length(ind.var)){
        AMIE.var.ind.mat.f[[i]] <- cbind(AMIE.var.ind[ind.var[[i]][,1]],AMIE.var.ind[ind.var[[i]][,2]])
    }
    
    ## Compute Conditional Effect
    if(inference==TRUE){
        CE.l <- list()
        for(i in 1:length(ind.var)){
            CE.tab <- matrix(NA, ncol=4, nrow=nrow(AMIE.var.ind.mat.f[[i]]))
            for(ce in 1:nrow(AMIE.var.ind.mat.f[[i]])){
                point <- (coefs.u[AME.var.ind.mat.f[[i]][ce,1]] - coefs.u[AME.var.ind.mat.f[[i]][ce,2]]) +
                    (coefs.u[AMIE.var.ind.mat.f[[i]][ce,1]] - coefs.u[AMIE.var.ind.mat.f[[i]][ce,2]]) 
                vari <- VarCondEffect(AME.var.ind.mat.f[[i]][ce,1], AME.var.ind.mat.f[[i]][ce,2],
                                      AMIE.var.ind.mat.f[[i]][ce,1], AMIE.var.ind.mat.f[[i]][ce,2], vcov.full=vcov.u)
                if(vari>0) std <- sqrt(vari) else std <- 0 
                CE.tab[ce,1:4] <- c(point, std, point - 1.96*std, point + 1.96*std)
            }
            rownames(CE.tab) <- level.name[[treat.ind]]
            colnames(CE.tab) <- c("ConditionalEffect", "sd", "2.5% CI", "97.5% CI")
            CE.tab <- round(CE.tab, digits=round)
            CE.l[[i]] <- CE.tab
        }
        names(CE.l) <- paste(var.name[cond.ind], "=", level.name[[cond.ind]],sep="")
    }else{
        CE.l <- list()
        for(i in 1:length(ind.var)){
            CE.tab <- c()
            for(ce in 1:nrow(AMIE.var.ind.mat.f[[i]])){
                point <- (coefs.u[AME.var.ind.mat.f[[i]][ce,1]] - coefs.u[AME.var.ind.mat.f[[i]][ce,2]]) +
                    (coefs.u[AMIE.var.ind.mat.f[[i]][ce,1]] - coefs.u[AMIE.var.ind.mat.f[[i]][ce,2]])                 
                CE.tab[ce] <- point
            }
            names(CE.tab) <- level.name[[treat.ind]]
            CE.tab <- round(CE.tab, digits=round)
            CE.l[[i]] <- CE.tab
        }
        names(CE.l) <- paste(var.name[cond.ind], "=", level.name[[cond.ind]],sep="")
    }
    
    
    ## ## base.ind
    ## ## if(missing(base.ind)==TRUE) base.ind <- 1
    
    ## AME.u0 <- AME[treat.ind][[1]] - AME[treat.ind][[1]][base.ind]
    ## AME.u <- matrix(rep(AME.u0, fac.level[cond.ind]), nrow=fac.level[cond.ind], byrow=TRUE)
    
    ## AMIE.mat <- matrix(unlist(AMIE2[INT.ind]), nrow=fac.level[Fac.ind[1]], ncol=fac.level[Fac.ind[2]])
    ## if(Norotate==FALSE){
    ##     AMIE.mat <- t(AMIE.mat)
    ## }else{
    ##         
    ## AMIE.u  <- AMIE.mat - AMIE.mat[,base.ind]
    
    ## CE <- AME.u + AMIE.u
    
    if(verbose==TRUE){
        cat(paste("\nTreatment Factor is ", var.name[treat.ind], " and ",
                  "Conditioning Factor is ", var.name[cond.ind],"\n",sep=""))
        print(CE.l)
    }
    
    output <- list("ConditionalEffects"=CE.l,
                   "treat.fac"=treat.fac,
                   "cond.fac"=cond.fac,
                   "treat.level"=treat.level,
                   "cond.level"=cond.level)
}




# CreateCoef --------------------------------------------------------------

CreateCoef <- function(coef.use, Index.use, Fac.level.use, Gorder){
    
    OneANOVA <- function(Ovec){
        OvecF<- c(Ovec,-sum(Ovec))
        return(OvecF)
    }
    TwoANOVA <- function(Tmatrix){
        TmatrixF <- cbind(Tmatrix,-apply(Tmatrix,1,sum))
        TmatrixF2 <- rbind(TmatrixF,-apply(TmatrixF,2,sum))
        return(TmatrixF2)
    }
    ThreeANOVA <- function(Tarray,array.dim){
        TarrayF <- array(NA,(array.dim+1))
        ## Fill everything
        for(k in 1:(array.dim[3])){
            if(is.numeric(Tarray[,,k])){
                Tarray.k <- matrix(Tarray[,,k],nrow=array.dim[1],ncol=array.dim[2])
            }else{
                Tarray.k <- Tarray[,,k]
            }
            TarrayF[,,k] <- TwoANOVA(Tarray.k)
        }
        ## Fill the last array
        coef.mat.sum <- matrix(0,nrow=nrow(TarrayF[,,1]),ncol=ncol(TarrayF[,,1]))
        for(k in 1:(array.dim[3])){
            coef.mat.sum <- coef.mat.sum + TarrayF[,,k]
        }        
        TarrayF[,,(array.dim[3]+1)] <- - coef.mat.sum
        return(TarrayF)
    }
    
    
    One.coef <- list()
    Two.coef <- list()
    Three.coef <- list()
    one <- two <- three <- 0
    for(z in 1:nrow(Index.use)){
        if(Index.use$order[z]==1){
            one <- one + 1
            One.coef[[one]] <- OneANOVA(coef.use[Index.use$start[z]:Index.use$end[z]])
        }
        if(Index.use$order[z]==2){
            two <- two + 1
            two.level <- Fac.level.use[z,][Fac.level.use[z,]!=0]
            two.mat <- matrix(coef.use[Index.use$start[z]:Index.use$end[z]],
                              nrow=two.level[1],ncol=two.level[2])
            two.mat <- matrix(two.mat,nrow=two.level[1],ncol=two.level[2])
            Two.coef[[two]] <- TwoANOVA(Tmatrix=two.mat)
        }
        if(Index.use$order[z]==3){
            three <- three + 1
            three.level <- Fac.level.use[z,][Fac.level.use[z,]!=0]
            three.array <- array(coef.use[Index.use$start[z]:Index.use$end[z]],
                                 dim=three.level)
            Three.coef[[three]] <- ThreeANOVA(three.array,three.level)
        }
    }
    if(Gorder ==1){
        full <- One.coef
    }else if(Gorder==2){
        full <- append(One.coef, Two.coef)
    }else if(Gorder==3){
        full <- append(One.coef, Two.coef)
        full <- append(full, Three.coef)
    }
    output <- list("One.coef"=One.coef,
                   "Two.coef"=Two.coef,
                   "Three.coef"=Three.coef,
                   "full"=full)
    return(output)
    
}

# CreateWeights -----------------------------------------------------------

CreateWeights <- function(formula,Gorder=3,dif=TRUE,facCons=FALSE,data,side,type.ind,fac.level,ord.fac,
                          indTwo=NULL, indThree=NULL){
    
    ## Make the difference in three-way interactions
    D3function <- function(alevel,blevel,clevel,type="A",ordABC=c(FALSE,FALSE,FALSE)){
        if(type=="A"){
            D3 <- D3Afunction(alevel=alevel,blevel=blevel,clevel=clevel,ord=ordABC[1])
        }else if(type=="B"){
            D3 <- D3Bfunction(alevel=alevel,blevel=blevel,clevel=clevel,ord=ordABC[2])
        }else if(type=="C"){
            D3 <- D3Cfunction(alevel=alevel,blevel=blevel,clevel=clevel,ord=ordABC[3])
        }
        return(D3)
    }
    
    D3Afunction <- function(alevel,blevel,clevel,ord=FALSE){
        D3a <- matrix(0,nrow=0,ncol=(alevel*blevel*clevel))  
        if(ord==FALSE){
            for(i in 1:(alevel-1)){
                for(j  in (i+1):alevel){
                    base.v <- rep(0,alevel)
                    base.v[i] <- -1
                    base.v[j] <- 1
                    base.mat <- diag((blevel*clevel)) %x% t(base.v)
                    D3a <- rbind(D3a,base.mat)
                }
            }
        }else if(ord==TRUE){
            for(i in 1:(alevel-1)){
                base.v <- rep(0,alevel)
                base.v[i] <- -1
                base.v[(i+1)] <- 1
                base.mat <- diag((blevel*clevel)) %x% t(base.v)
                D3a <- rbind(D3a,base.mat)
            }        
        }
        return(D3a)
    }
    
    ## (B): Make the difference in two-way interactions 
    D3Bfunction <- function(alevel,blevel,clevel,ord=FALSE){    
        if(ord==FALSE){
            D3b <- matrix(0,nrow=0,ncol=(alevel*blevel*clevel))
            for(i in 1:(blevel-1)){     
                left.temp <- matrix(0,ncol=((i-1)*(alevel)),nrow=(alevel))
                left.main <- -diag((alevel))
                left.final <- cbind(left.temp,left.main)
                for(j  in (i+1):blevel){
                    middle.temp <-  matrix(0,ncol=((j-i-1)*(alevel)),nrow=(alevel))
                    right.main <- diag((alevel))
                    right.temp <- matrix(0,ncol=((blevel-j)*(alevel)),nrow=(alevel))
                    right.final <- cbind(middle.temp,right.main,right.temp)
                    D2b <- cbind(left.final,right.final) ## The base
                    D3base <- diag(clevel) %x% D2b
                    D3b <- rbind(D3b, D3base)
                }
            }
        }else if (ord==TRUE){
            D3b <- matrix(0,nrow=0,ncol=(alevel*blevel*clevel))
            for(i in 1:(blevel-1)){     
                left <- matrix(0,ncol=((i-1)*(alevel)),nrow=(alevel))
                Main <- cbind(-diag((alevel)),diag((alevel)))
                right <- matrix(0,ncol=((blevel-(i+1))*(alevel)),nrow=(alevel))
                D2b <- cbind(left,Main,right)
                D3base <- diag(clevel) %x% D2b
                D3b <- rbind(D3b, D3base)
            }
        }
        return(D3b)
        ## psy_B(2,3), psy_B(2,4), psy_B(2,5) 
    }
    
    D3Cfunction <- function(alevel,blevel,clevel,ord=FALSE){
        if(ord==FALSE){
            D3c <- matrix(0,nrow=0,ncol=(alevel*blevel*clevel))
            for(i in 1:(clevel-1)){     
                left.temp <- matrix(0,ncol=((i-1)*(alevel*blevel)),nrow=(alevel*blevel))
                left.main <- -diag((alevel*blevel))
                left.final <- cbind(left.temp,left.main)
                for(j  in (i+1):clevel){
                    middle.temp <-  matrix(0,ncol=((j-i-1)*(alevel*blevel)),nrow=(alevel*blevel))
                    right.main <- diag((alevel*blevel))
                    right.temp <- matrix(0,ncol=((clevel-j)*(alevel*blevel)),
                                         nrow=(alevel*blevel))
                    right.final <- cbind(middle.temp,right.main,right.temp)
                    D3c <- rbind(D3c,cbind(left.final,right.final))
                }
            }
        }else if (ord==TRUE){
            Left <- cbind(-diag((alevel*blevel)*(clevel-1)),
                          matrix(0,ncol=(alevel*blevel),nrow=(alevel*blevel)*(clevel-1)))
            Right <- cbind(matrix(0,ncol=(alevel*blevel),nrow=(alevel*blevel)*(clevel-1)),
                           diag((alevel*blevel)*(clevel-1)))
            D3c <- Left + Right
        }
        return(D3c)
    }
    
    
    ## From this version, I use D2 rather than D3 
    ## Make the difference in two-way interactions 
    D2function <- function(alevel,blevel,type,ordAB){
        if(type=="A"){
            D2 <- D2Afunction(alevel=alevel,blevel=blevel,ord=ordAB[1])
        }else if(type=="B"){        
            D2 <- D2Bfunction(alevel=alevel,blevel=blevel,ord=ordAB[2])
        }
        return(D2)
        ## D2 matrix
        ## nrow = differences between two-way interactions fixing the sublevel
        ## ncol = alevel*blevel (coefficients for two-ways)
    }
    
    ## (A): Make the difference in two-way interactions 
    D2Afunction <- function(alevel,blevel,ord=FALSE){
        D2a <- matrix(0,nrow=0,ncol=(alevel*blevel))
        if(ord==FALSE){
            for(i in 1:(alevel-1)){
                for(j  in (i+1):alevel){
                    base.v <- rep(0,alevel)
                    base.v[i] <- -1
                    base.v[j] <- 1
                    base.mat <- diag(blevel) %x% t(base.v)
                    D2a <- rbind(D2a,base.mat)
                }
            }
        }else if(ord==TRUE){
            for(i in 1:(alevel-1)){
                base.v <- rep(0,alevel)
                base.v[i] <- -1
                base.v[(i+1)] <- 1
                base.mat <- diag(blevel) %x% t(base.v)
                D2a <- rbind(D2a,base.mat)
            }        
        }
        return(D2a)
    }
    
    ## (B): Make the difference in two-way interactions 
    D2Bfunction <- function(alevel,blevel,ord=FALSE){    
        if(ord==FALSE){
            D2b <- matrix(0,nrow=0,ncol=(alevel*blevel))
            for(i in 1:(blevel-1)){     
                left.temp <- matrix(0,ncol=((i-1)*(alevel)),nrow=(alevel))
                left.main <- -diag((alevel))
                left.final <- cbind(left.temp,left.main)
                for(j  in (i+1):blevel){
                    middle.temp <-  matrix(0,ncol=((j-i-1)*(alevel)),nrow=(alevel))
                    right.main <- diag((alevel))
                    right.temp <- matrix(0,ncol=((blevel-j)*(alevel)),nrow=(alevel))
                    right.final <- cbind(middle.temp,right.main,right.temp)
                    D2b <- rbind(D2b,cbind(left.final,right.final))
                }
            }
        }else if (ord==TRUE){
            Left <- cbind(-diag((alevel)*(blevel-1)),matrix(0,ncol=alevel,nrow=(alevel)*(blevel-1)))
            Right <- cbind(matrix(0,ncol=alevel,nrow=(alevel)*(blevel-1)), diag((alevel)*(blevel-1)))
            D2b <- Left + Right
        }
        return(D2b)
        ## psy_B(2,3), psy_B(2,4), psy_B(2,5) 
    }
    
    ## (Main): Making the difference between all levels
    D1function <- function(nlevel,ord=FALSE){
        if(ord==FALSE){
            qp = (nlevel)*(nlevel-1)/2
            D1.mat <- matrix(0,ncol=(nlevel),nrow=0)    
            if (qp==1) D1.mat <- matrix(c(-1,1),ncol=2,nrow=1)
            if (qp>1){
                for(i in 1 : (nlevel-1)){
                    w.diag <- diag((nlevel-i))
                    left.v <- c(rep(0,(i-1)),-1)
                    left.mat <- matrix(rep(left.v,(nlevel-i)),nrow=(nlevel-i),byrow=TRUE)
                    w.mat <- cbind(left.mat,w.diag)
                    D1.mat <- rbind(D1.mat,w.mat)
                }
            }
        }else if(ord==TRUE){
            if(nlevel==2){
                D1.mat <- matrix(c(-1,1),ncol=2,nrow=1)
            }else{
                D1.mat <- cbind(0,diag(nlevel-1)) + cbind(-diag(nlevel-1),0)
            }
            
        }    
        return(D1.mat)
        ## nrow = differences between main effects
        ## ncol = mainlevel (coefficients for the main factor)
    }
    
    if(dif==TRUE){
        data1 <- data[side==1,]
        data2 <- data[side==0,]
        model.frame  <- model.frame(formula,data=data1)
        contr <- rep(list("contr.sum"), ncol(model.frame) - 1)
        names(contr) <- colnames(model.frame)[-1]
        y <- model.frame[,1]
        x1 <- model.matrix(formula, data=data1,contrast=contr)[,-1]
        x2 <- model.matrix(formula, data=data2,contrast=contr)[,-1]
        X <- x1-x2
    }else if (dif==FALSE){
        model.frame  <- model.frame(formula,data=data)
        contr <- rep(list("contr.sum"), ncol(model.frame) - 1)
        names(contr) <- colnames(model.frame)[-1]
        y <- model.frame[,1]
        X <- model.matrix(formula, data=data,contrast=contr)[,-1]
    }
    lm.fit <- lm(y ~ X)
    coef.use <- coef(lm.fit)[-1]
    
    levelIndex <- CreatelevelIndex(fac.level=(fac.level-1),ord.fac=ord.fac,Gorder=Gorder,
                                   indTwo=indTwo, indThree=indThree)
    use.ind <-  (levelIndex$plus==1)*(levelIndex$dif==0)
    Index.use <- levelIndex[use.ind==1,]
    Fac.index <- levelIndex[,regexpr("Fac",colnames(levelIndex))>0]
    Fac.Ind.use <- Fac.index[use.ind==1,]
    Fac.level.useC <- (rep(1,nrow(Fac.Ind.use)) %x% t((fac.level-1))) * Fac.Ind.use
    Fac.level.use <- (rep(1,nrow(Fac.Ind.use)) %x% t((fac.level))) * Fac.Ind.use
    Index.use$start <- c(1,cumsum(Index.use$length)[-nrow(Index.use)]+1)
    Index.use$end <-  cumsum(Index.use$length)
    
    if(Gorder==2){
        Index.use$seq.order <- c(seq(1:sum(Index.use$order==1)),
                                 seq(1:sum(Index.use$order==2)))
    }else if(Gorder==3){
        Index.use$seq.order <- c(seq(1:sum(Index.use$order==1)),
                                 seq(1:sum(Index.use$order==2)),
                                 seq(1:sum(Index.use$order==3)))
    }
    
    ## Need to change?
    Coef.list <- CreateCoef(coef.use=coef.use, 
                            Index.use=Index.use,
                            Fac.level.use=Fac.level.useC,
                            Gorder=Gorder)
    
    if(facCons==TRUE){
        n.fac <- length(fac.level)
        MainDif <- list()
        IntDif <- list()
        ThreeDif <- list()
        for(z in 1:n.fac){
            Maincoef <- Coef.list$One.coef[[z]]
            MainD1 <- D1function(nlevel=fac.level[z],ord=ord.fac[z])
            MainDif[[z]] <- t(MainD1 %*% Maincoef)
        }
        
        ind.use.Two <- (Index.use$order==2)
        Index.useTwo <- Index.use[ind.use.Two==1,]
        Fac.useTwo <- Fac.Ind.use[ind.use.Two==1,]
        Fac.level.useTwo <- Fac.level.use[ind.use.Two==1,]
        ## IntDif <- matrix(0,nrow=0,ncol=qp)
        for(z in 1:nrow(Index.useTwo)){
            level.use <- as.numeric(Fac.level.useTwo[z,Fac.useTwo[z, ]==1])
            ord.use <- ord.fac[Fac.useTwo[z, ]==1]
            coef.int <- c(Coef.list$Two.coef[[z]])
            D2A  <- D2function(alevel=level.use[1],blevel=level.use[2],
                               type="A",ordAB=ord.use)
            D2B  <- D2function(alevel=level.use[1],blevel=level.use[2],
                               type="B",ordAB=ord.use)
            IntDif.A <- D2A %*% coef.int
            IntDif.B <- D2B %*% coef.int
            IntDif[[z]] <- c(IntDif.A, IntDif.B)
            
            if(Gorder>2){
                ## Three-ways   
                ind.use.Three <- (Index.use$order==3)
                Index.useThree <- Index.use[ind.use.Three==1,]
                Fac.useThree <- Fac.Ind.use[ind.use.Three==1,]
                Fac.level.useThree <- Fac.level.use[ind.use.Three==1,]
                ThreeDif <- list()
                
                for(z in 1:nrow(Index.useThree)){
                    level.use <- as.numeric(Fac.level.useThree[z,Fac.useThree[z, ]==1])
                    ord.use <- ord.fac[Fac.useThree[z, ]==1]                    
                    D3A  <- D3function(alevel=level.use[1],blevel=level.use[2],clevel=level.use[3],
                                       type="A",ordABC=ord.use)
                    D3B  <- D3function(alevel=level.use[1],blevel=level.use[2],clevel=level.use[3],
                                       type="B",ordABC=ord.use)
                    D3C  <- D3function(alevel=level.use[1],blevel=level.use[2],clevel=level.use[3],
                                       type="C",ordABC=ord.use)
                    coef.three <- c(Coef.list$Three.coef[[z]])
                    ThreeDif.A <- D3A %*% coef.three
                    ThreeDif.B <- D3B %*% coef.three
                    ThreeDif.C <- D3C %*% coef.three
                    ThreeDif[[z]] <- c(ThreeDif.A, ThreeDif.B, ThreeDif.C)
                }
            }            
        }
        
        Main.max <- unlist(lapply(1:length(MainDif),function(x) 1/max(abs(MainDif[[x]]))))
        Int.max <- unlist(lapply(1:length(IntDif),function(x) 1/max(abs(IntDif[[x]]))))        
        if(Gorder > 2){
            Three.max <- unlist(lapply(1:length(ThreeDif),function(x) 1/max(abs(ThreeDif[[x]]))))
            weight <- c(Main.max,Int.max, Three.max)
        }else if(Gorder==2){
            weight <- c(Main.max,Int.max)
        }
    }
    else if(facCons==FALSE){
        ## Main Effects 
        Maincoef <- Coef.list$One.coef[[type.ind]]
        MainD1 <- D1function(nlevel=fac.level[type.ind],ord=ord.fac[type.ind])
        MainDif <- t(MainD1 %*% Maincoef)
        qp <- ncol(MainDif)
        
        ## Two-ways
        onlyOne <- sum(Fac.Ind.use[type.ind])==1
        yesTwo <- sum((Index.use$order==2)*(Fac.Ind.use[type.ind]==1))>0
        yesThree <- sum((Index.use$order==3)*(Fac.Ind.use[type.ind]==1))>0
        if(onlyOne==TRUE){
            ## No Interaction
            Dif <- abs(MainDif)
            ## print(Dif)
        }else if(yesTwo==TRUE){
            ind.use.Two <- (Index.use$order==2)*(Fac.Ind.use[type.ind]==1)
            Index.useTwo <- Index.use[ind.use.Two==1,]
            Fac.useTwo <- Fac.Ind.use[ind.use.Two==1,]
            Fac.level.useTwo <- Fac.level.use[ind.use.Two==1,]
            IntDif <- matrix(0,nrow=0,ncol=qp)
            for(z in 1:nrow(Index.useTwo)){
                level.use <- as.numeric(Fac.level.useTwo[z,Fac.useTwo[z, ]==1])
                ord.use <- ord.fac[Fac.useTwo[z, ]==1]
                if(min(which(Fac.useTwo[z,]==1))==type.ind){
                    type.d2 <- "A"}else{ type.d2 <- "B" }
                D2  <- D2function(alevel=level.use[1],blevel=level.use[2],
                                  type=type.d2,ordAB=ord.use)
                coef.int <- c(Coef.list$Two.coef[[Index.useTwo[z,"seq.order"]]])
                IntDif.t <- D2 %*% coef.int
                IntDif.mat <- matrix(IntDif.t,ncol=qp)
                IntDif <- rbind(IntDif, IntDif.mat)
            }
            if(yesThree==FALSE){                
                Dif <- abs(rbind(MainDif,IntDif))
            }else if(yesThree==TRUE){
                ## Three-ways   
                ind.use.Three <- (Index.use$order==3)*(Fac.Ind.use[type.ind]==1)
                Index.useThree <- Index.use[ind.use.Three==1,]
                Fac.useThree <- Fac.Ind.use[ind.use.Three==1,]
                Fac.level.useThree <- Fac.level.use[ind.use.Three==1,]
                ThreeDif <- matrix(0,nrow=0,ncol=qp)
                
                for(z in 1:nrow(Index.useThree)){
                    level.use <- as.numeric(Fac.level.useThree[z,Fac.useThree[z, ]==1])
                    ord.use <- ord.fac[Fac.useThree[z, ]==1]
                    if(min(which(Fac.useThree[z,]==1))==type.ind){
                        type.d3 <- "A"
                    }else if (max(which(Fac.useThree[z,]==1))==type.ind){
                        type.d3 <- "C"
                    }else{ type.d3 <- "B"}
                    D3  <- D3function(alevel=level.use[1],blevel=level.use[2],clevel=level.use[3],
                                      type=type.d3,ordABC=ord.use)
                    coef.three <- c(Coef.list$Three.coef[[Index.useThree[z,"seq.order"]]])
                    ThreeDif.t <- D3 %*% coef.three
                    ThreeDif.mat <- matrix(ThreeDif.t,ncol=qp)
                    ThreeDif <- rbind(ThreeDif, ThreeDif.mat)
                }            
                Dif <- abs(rbind(MainDif,IntDif,ThreeDif))
            }
        }
        
        ## if(Gorder == 2){
        ##     Dif <- abs(rbind(MainDif,IntDif))
        ## }else if(Gorder > 2){
        ##     Dif <- abs(rbind(MainDif,IntDif,ThreeDif))
        ## }
        weight <- apply(Dif,2,function(x) 1/max(abs(x)))
    }
    
    levelIndex <- CreatelevelIndex(fac.level=(fac.level),ord.fac=ord.fac, Gorder=Gorder,
                                   indTwo=indTwo, indThree=indThree)
    use.ind <-  (levelIndex$plus==1)*(levelIndex$dif==1) * (levelIndex$order==1)
    Index.use <- levelIndex[use.ind==1,]
    temp.w <- 1/(sqrt(fac.level[type.ind])*(fac.level[type.ind]+1))
    weight.u <- rep(temp.w, times=Index.use$length[type.ind])
    
    ## Combine two weights. 
    output <- list("weight.ols"=weight, "weight.fac"=weight.u)
    return(output)
}


# CreatelevelIndex --------------------------------------------------------

CreatelevelIndex <- function(fac.level,ord.fac, indTwo=NULL, indThree=NULL, Gorder){
    if(missing(Gorder)){
        Gorder <- 3 
    }
    n.fac <- length(fac.level)
    
    dif.level <- rep(0,times=n.fac)
    for(i in 1:n.fac){
        if(ord.fac[i]==TRUE){
            dif.level[i] <- fac.level[i] - 1
        }else{
            dif.level[i] <- fac.level[i]*(fac.level[i] - 1)/2
        }
    }
    
    fac.name <- paste("Fac",seq(1:n.fac),sep=".")
    
    ## Setup
    if(is.null(indTwo)==TRUE){
        indTwo <- combn(seq(1:n.fac),2)
    }
    if(is.null(indThree)==TRUE & Gorder==3){
        indThree <- combn(seq(1:n.fac),3)
    }
    
    ## %%%% IndTwo can be imported. 
    
    ## For one-way
    Index.one <- matrix(0,nrow=0,ncol=(5 + n.fac))
    BaseOne <- as.data.frame(cbind(rep(c(1,0),2),rep(c(0,1),each=2)))
    BaseOne <- cbind(rep(1,4),BaseOne)
    for(one in 1:n.fac){
        FacOne  <- matrix(0,ncol=n.fac,nrow=4)
        FacOne[,one] <- rep(c(1,2),each=2)
        length <- c(rep(fac.level[one],each=2),rep(dif.level[one],each=2))
        Index.one <- rbind(Index.one, cbind(length,BaseOne,FacOne))
    }
    colnames(Index.one) <- c("length","order","plus","dif",fac.name)
    ## Index.one
    
    ## For two-way
    if(Gorder>=2){
        Index.two <- matrix(0,nrow=0,ncol=(5 + n.fac))
        BaseTwo <- as.data.frame(cbind(rep(c(1,0),3),c(rep(0,2),rep(1,4))))
        BaseTwo <- cbind(rep(2,6),BaseTwo)
        for(z in 1:ncol(indTwo)){
            FacTwo  <- matrix(0,ncol=n.fac,nrow=6)
            FacTwo[,indTwo[,z]] <- cbind(rep(c(1,2,1),each=2),rep(c(1,1,2),each=2))
            fac.int.level <- fac.level[indTwo[1,z]]*fac.level[indTwo[2,z]]
            fac.int1 <- dif.level[indTwo[1,z]]*fac.level[indTwo[2,z]]
            fac.int2 <- fac.level[indTwo[1,z]]*dif.level[indTwo[2,z]]
            length <- c(rep(fac.int.level,each=2),rep(fac.int1,each=2),rep(fac.int2,each=2))
            Index.two <- rbind(Index.two, cbind(length,BaseTwo,FacTwo))
        }
        colnames(Index.two) <- c("length","order","plus","dif",fac.name)
    }
    ## Index.two
    
    if(Gorder == 3){
        ## For three-way
        Index.three <- matrix(0,nrow=0,ncol=(5 + n.fac))
        BaseThree <- as.data.frame(cbind(rep(c(1,0),4),c(rep(0,2),rep(1,6))))
        BaseThree <- cbind(rep(3,8),BaseThree)
        for(z in 1:ncol(indThree)){
            FacThree  <- matrix(0,ncol=n.fac,nrow=8)
            FacThree[,indThree[,z]] <- cbind(rep(c(1,2,1,1),each=2),
                                             rep(c(1,1,2,1),each=2),
                                             rep(c(1,1,1,2),each=2))
            fac.int.level <- fac.level[indThree[1,z]]*fac.level[indThree[2,z]]*fac.level[indThree[3,z]]
            fac.int1 <- dif.level[indThree[1,z]]*fac.level[indThree[2,z]]*fac.level[indThree[3,z]]
            fac.int2 <- fac.level[indThree[1,z]]*dif.level[indThree[2,z]]*fac.level[indThree[3,z]]
            fac.int3 <- fac.level[indThree[1,z]]*fac.level[indThree[2,z]]*dif.level[indThree[3,z]]
            length <- c(rep(fac.int.level,each=2),rep(fac.int1,each=2),
                        rep(fac.int2,each=2), rep(fac.int3,each=2))
            Index.three <- rbind(Index.three, cbind(length,BaseThree,FacThree))
        }
        colnames(Index.three) <- c("length","order","plus","dif",fac.name)
    }
    
    ## Add Start Index
    if(Gorder==1){
        Index.Mat <- Index.one
    }else if(Gorder ==2){
        Index.Mat <- rbind(Index.one, Index.two)
    }else if(Gorder ==3){
        Index.Mat <- rbind(Index.one, Index.two,Index.three)
    }
    start <- c(1,cumsum(Index.Mat$length)[-nrow(Index.Mat)]+1)
    
    Index.Mat <- cbind(start,Index.Mat)
    
    return(Index.Mat)
}




# FullVCOV ----------------------------------------------------------------

IndexOne <- function(index, n.level){
    index.u <- list()
    base.u <- rep(0, n.level)
    if(index <=n.level - 1){
        index.u[[1]] <- index   ## just an index in the trancated scale 
        base.u[index] <- 1      ## matrix form 
        index.u[[2]] <- base.u  ## matrix 
        index.u[[3]] <- 1  ## plus
    }else{
        index.u[[1]] <- seq(1:(n.level-1))
        base.u[-index] <- -1
        index.u[[2]] <- base.u  ## matrix 
        index.u[[3]] <- -1  ## plus
    }
    return(index.u)
}

IndexTwo <- function(index1, index2, n.level1, n.level2){
    index.u <- list()
    base.u <- matrix(0, nrow=n.level1, ncol=n.level2)
    if(index1 <= n.level1 - 1 & index2 <= n.level2 - 1){
        index.u[[1]] <- index1 + (index2 - 1)*(n.level1 - 1)   ## just an index in the trancated scale 
        base.u[index1, index2] <- 1      ## matrix form 
        index.u[[2]] <- base.u  ## matrix 
        index.u[[3]] <- 1  ## plus
    }else if(index1 == n.level1 & index2 <= n.level2 -1){
        index.u[[1]] <- seq(from=(1 + (index2 - 1)*(n.level1 - 1)), to=(index2*(n.level1 - 1)))
        base.u[1:(n.level1-1), index2] <- -1
        index.u[[2]] <- base.u  ## matrix 
        index.u[[3]] <- -1  ## plus
    }else if(index1 <= (n.level1 -1) & index2 == n.level2){
        index.u[[1]] <- seq(from=index1, to=(index1+(n.level1 - 1)*(n.level2-2)), by=(n.level1-1))
        base.u[index1, 1:(n.level2-1)] <- -1
        index.u[[2]] <- base.u  ## matrix 
        index.u[[3]] <- -1  ## plus
    }else if(index1 == n.level1 & index2 == n.level2){
        index.u[[1]] <- seq(from=1, to=(n.level1 - 1)*(n.level2-1))
        base.u[1:(n.level1-1), 1:(n.level2-1)] <- 1
        index.u[[2]] <- base.u  ## matrix 
        index.u[[3]] <- 1  ## plus
    }
    return(index.u)
}

IndexThree <- function(index1, index2, index3, n.level1, n.level2, n.level3){
    index.u <- list()
    base.u <- array(0, dim=c(n.level1,n.level2, n.level3))
    if(index3 <= n.level3 - 1){
        if(index1 <= n.level1 - 1 & index2 <= n.level2 - 1){
            index.u[[1]] <- index1 + (index2 - 1)*(n.level1 - 1) +
                (n.level1-1)*(n.level2-1)*(index3-1) ## just an index in the trancated scale 
            base.u[index1, index2, index3] <- 1      ## matrix form 
            index.u[[2]] <- base.u  ## matrix 
            index.u[[3]] <- 1  ## plus
        }else if(index1 == n.level1 & index2 <= n.level2 -1){
            index.u[[1]] <- seq(from=(1 + (index2 - 1)*(n.level1 - 1)), to=(index2*(n.level1 - 1))) +
                (n.level1-1)*(n.level2-1)*(index3-1)
            base.u[1:(n.level1-1), index2, index3] <- -1
            index.u[[2]] <- base.u  ## matrix 
            index.u[[3]] <- -1  ## plus
        }else if(index1 <= (n.level1 -1) & index2 == n.level2){
            index.u[[1]] <- seq(from=index1, to=(index1+(n.level1 - 1)*(n.level2-2)), by=(n.level1-1)) +
                (n.level1-1)*(n.level2-1)*(index3-1)
            base.u[index1, 1:(n.level2-1), index3] <- -1
            index.u[[2]] <- base.u  ## matrix 
            index.u[[3]] <- -1  ## plus
        }else if(index1 == n.level1 & index2 == n.level2){
            index.u[[1]] <- seq(from=1, to=(n.level1 - 1)*(n.level2-1)) +
                (n.level1-1)*(n.level2-1)*(index3-1)
            base.u[1:(n.level1-1), 1:(n.level2-1), index3] <- 1
            index.u[[2]] <- base.u  ## matrix 
            index.u[[3]] <- 1  ## plus
        }
    }else if(index3==n.level3){
        if(index1 <= n.level1 - 1 & index2 <= n.level2 - 1){
            index.u[[1]] <- seq(from=(index1 + (index2 - 1)*(n.level1 - 1)),
                                to=((index1 + (index2 - 1)*(n.level1 - 1)) +
                                        (n.level1-1)*(n.level2-1)*(n.level3-2)),
                                by=(n.level1-1)*(n.level2-1))  ## just an index in the trancated scale 
            base.u[index1, index2, seq(1:(n.level3-1))] <- -1      ## matrix form 
            index.u[[2]] <- base.u  ## matrix 
            index.u[[3]] <- -1  ## plus
        }else if(index1 == n.level1 & index2 <= n.level2 -1){
            b.i <- seq(from=(1 + (index2 - 1)*(n.level1 - 1)), to=(index2*(n.level1 - 1)))
            b.i.l <- c()
            for(z in 1:(n.level3-1)){
                b.i.l <- c(b.i.l, (n.level1-1)*(n.level2-1)*(z-1) + b.i)
            }
            index.u[[1]] <- b.i.l
            base.u[1:(n.level1-1), index2, seq(1:(n.level3-1))] <- 1
            index.u[[2]] <- base.u  ## matrix 
            index.u[[3]] <- 1  ## plus
        }else if(index1 <= (n.level1 -1) & index2 == n.level2){
            b.i <- seq(from=index1, to=(index1+(n.level1 - 1)*(n.level2-2)), by=(n.level1-1))
            b.i.l <- c()
            for(z in 1:(n.level3-1)){
                b.i.l <- c(b.i.l, (n.level1-1)*(n.level2-1)*(z-1) + b.i)
            }
            index.u[[1]] <- b.i.l
            base.u[index1, 1:(n.level2-1), seq(1:(n.level3-1))] <- 1
            index.u[[2]] <- base.u  ## matrix 
            index.u[[3]] <- -1  ## plus
        }else if(index1 == n.level1 & index2 == n.level2){
            index.u[[1]] <- seq(from=1, to=(n.level1 - 1)*(n.level2-1)*(n.level3-1))
            base.u[1:(n.level1-1), 1:(n.level2-1), seq(1:(n.level3-1))] <- -1
            index.u[[2]] <- base.u  ## matrix 
            index.u[[3]] <- -1  ## plus
        }
    }
    return(index.u)
}

IndConvert <- function(index.level, index.fac, term.index,
                       fac.level, Index.use){
    nway <- length(index.level)
    fac.use <- fac.level[index.fac]
    if(term.index==1) base.n <- 0 else base.n <- Index.use$end[(term.index-1)]
    if(nway==1) base <- IndexOne(index=index.level, n.level=fac.use)
    if(nway==2) base <- IndexTwo(index1=index.level[1], index2=index.level[2],
                                 n.level1=fac.use[1], n.level2=fac.use[2])
    if(nway==3) base <- IndexThree(index1=index.level[1], index2=index.level[2], index3=index.level[3],
                                   n.level1=fac.use[1], n.level2=fac.use[2], n.level3=fac.use[3])
    base.ind <- base[[1]]
    sign <- base[[3]]
    ## final ind
    index.u <- base.n + base.ind
    out <- list("index.u"=index.u, "sign"=sign)
    return(out)
}

ind2Var <- function(out, vcov){
    ## index.u from IndConvert
    index.u <- out$index.u
    if(length(index.u)==1){
        Var <- vcov[index.u, index.u]
    }else{
        ## Linear Combination
        Var <- sum(vcov[index.u, index.u])
    }
    return(Var)
}

ind2Cov <- function(out1, out2, vcov){
    ## index.u1 and index.u2 from IndConvert
    index.u1 <- out1$index.u
    index.u2 <- out2$index.u
    sign1 <- matrix(rep(out1$sign, length(index.u1)),nrow=1)
    sign2 <- matrix(rep(out2$sign, length(index.u2)),nrow=1)        
    base.COV <- vcov[index.u1, index.u2]
    Cov <- sign1 %*% base.COV %*% t(sign2)
    return(Cov)
}
## Use these to construct FULL Variance Matrix
## I want to fill in based on rows.

VarEffect <- function(ind1, ind2, vcov.full){    
    effect.var <- vcov.full[ind1,ind1] + vcov.full[ind2,ind2] - 2*vcov.full[ind1,ind2]
    return(effect.var)
}

VarCondEffect <- function(ame.ind1, ame.ind2,
                          amie.ind1, amie.ind2,
                          vcov.full){
    if(ame.ind1== ame.ind2 & amie.ind1==amie.ind2){
        effect.var <- 0
    }else{
        effect.var <- vcov.full[ame.ind1,ame.ind1] + vcov.full[ame.ind2,ame.ind2] +
            vcov.full[amie.ind1,amie.ind1] + vcov.full[amie.ind2,amie.ind2] -
            2*vcov.full[ame.ind1,ame.ind2] + 2*vcov.full[ame.ind1,amie.ind1] -
            2*vcov.full[ame.ind1,amie.ind2] - 2*vcov.full[ame.ind2,amie.ind1] +
            2*vcov.full[ame.ind2,amie.ind2] - 2*vcov.full[amie.ind1,amie.ind2]
    }
    return(effect.var)
}

FullVCOV <- function(vcov, fac.level, ord.fac, indTwo, indThree, Gorder=Gorder){
    ## Now I need to conect Effects to Index.
    levelIndexS <- CreatelevelIndex(fac.level=(fac.level-1),ord.fac=ord.fac, 
                                    Gorder=Gorder, indTwo=indTwo, indThree=indThree)
    use.indS <-  (levelIndexS$plus==1)*(levelIndexS$dif==0)
    Index.useS <- levelIndexS[use.indS==1,]
    levelIndex <- CreatelevelIndex(fac.level=fac.level,ord.fac=ord.fac, 
                                   Gorder=Gorder, indTwo=indTwo, indThree=indThree)
    Index.useS$end <-  cumsum(Index.useS$length)
    use.ind <-  (levelIndex$plus==1)*(levelIndex$dif==0)
    Index.use <- levelIndex[use.ind==1,]    
    Fac.index <- levelIndex[,regexpr("Fac",colnames(levelIndex))>0]
    Fac.Ind.use <- Fac.index[use.ind==1,]
    
    out <- list()
    sign <- c()
    for(z in 1:nrow(Index.use)){
        if(Index.use$order[z]==1){
            out0 <- list()
            sign0 <- c()
            index.fac <- which(Fac.Ind.use[z,]==1)
            for(i in 1:Index.use$length[z]){
                out0[[i]] <- IndConvert(index.level=i, index.fac=index.fac, term.index=z,
                                        fac.level=fac.level, Index.use=Index.useS)
            }
            out <- c(out, out0)
        }
        if(Index.use$order[z]==2){
            out0 <- list()
            sign0 <- c()
            index.fac <- which(Fac.Ind.use[z,]==1)
            exp2 <- as.matrix(expand.grid(seq(1:fac.level[index.fac[1]]), seq(1:fac.level[index.fac[2]])))
            for(i in 1:Index.use$length[z]){
                out0[[i]] <- IndConvert(index.level=c(exp2[i,1:2]), index.fac=c(index.fac), term.index=z,
                                        fac.level=fac.level, Index.use=Index.useS)                
            }
            out <- c(out, out0)
        }             
        if(Index.use$order[z]==3){
            out0 <- list()
            index.fac <- which(Fac.Ind.use[z,]==1)
            exp3 <- as.matrix(expand.grid(seq(1:fac.level[index.fac[1]]),
                                          seq(1:fac.level[index.fac[2]]),
                                          seq(1:fac.level[index.fac[3]])))
            for(i in 1:Index.use$length[z]){
                out0[[i]] <- IndConvert(index.level=c(exp3[i,1:3]), index.fac=c(index.fac), term.index=z,
                                        fac.level=fac.level, Index.use=Index.useS)
                
            }
            out <- c(out, out0)               
        }
    }
    
    vcov.f <- matrix(0, ncol=sum(Index.use$length), nrow=sum(Index.use$length))
    vcov.f[1,1] <- ind2Var(out=out[[1]], vcov=vcov)
    for(i in 2:length(out)){           
        main.out <- out[[i]]
        for(j in 1:(i-1)){
            vcov.f[i,j] <- ind2Cov(out1=main.out, out2=out[[j]], vcov=vcov)
        }
        vcov.f[i,i] <- ind2Var(out=main.out, vcov=vcov)
    }
    vcov.f2 <- vcov.f
    diag(vcov.f2) <- 0
    vcov.out <- vcov.f + t(vcov.f2)
    
    return(vcov.out)
}


# NoRegularization --------------------------------------------------------

NoRegularization <- function(y, X.no, base.name, fac.level, ord.fac, Gorder, indTwo, indThree,
                             cluster){
    
    lm.fit.no <- lm(y ~ X.no)
    intercept <- coef(lm.fit.no)[1]
    coef.no <- coef(lm.fit.no)[-1]
    levelIndex <- CreatelevelIndex(fac.level=(fac.level-1),ord.fac=ord.fac,Gorder=Gorder,
                                   indTwo=indTwo, indThree=indThree)
    use.ind <-  (levelIndex$plus==1)*(levelIndex$dif==0)
    Index.use <- levelIndex[use.ind==1,]
    Fac.index <- levelIndex[,regexpr("Fac",colnames(levelIndex))>0]
    Fac.Ind.use <- Fac.index[use.ind==1,]
    Fac.level.useC <- (rep(1,nrow(Fac.Ind.use)) %x% t((fac.level-1))) * Fac.Ind.use
    Fac.level.use <- (rep(1,nrow(Fac.Ind.use)) %x% t((fac.level))) * Fac.Ind.use
    Index.use$start <- c(1,cumsum(Index.use$length)[-nrow(Index.use)]+1)
    Index.use$end <-  cumsum(Index.use$length)
    
    if(Gorder==1){
        Index.use$seq.order <- seq(1:sum(Index.use$order==1))
    }else if(Gorder==2){
        Index.use$seq.order <- c(seq(1:sum(Index.use$order==1)),
                                 seq(1:sum(Index.use$order==2)))
    }else if(Gorder==3){
        Index.use$seq.order <- c(seq(1:sum(Index.use$order==1)),
                                 seq(1:sum(Index.use$order==2)),
                                 seq(1:sum(Index.use$order==3)))
    }
    
    coefs.f <- CreateCoef(coef.use=coef.no, 
                          Index.use=Index.use,
                          Fac.level.use=Fac.level.useC,
                          Gorder=Gorder)$full
    
    if(is.null(cluster)==FALSE){
        vcov.c <- cluster_se_glm(model=lm.fit.no, cluster=cluster)
        vcov.u <- vcov.c[-1, -1]
    }else{
        vcov.u <- vcov(lm.fit.no)[-1,-1]
    }
    
    
    vcov.f <- FullVCOV(vcov=vcov.u, fac.level=fac.level, ord.fac=ord.fac,
                       indTwo=indTwo, indThree=indThree, Gorder=Gorder)
    
    ## Name
    levelIndex2 <- CreatelevelIndex(fac.level=fac.level,ord.fac=ord.fac,Gorder=Gorder,
                                    indTwo=indTwo, indThree=indThree)
    use.ind2 <-  (levelIndex2$plus==1)*(levelIndex2$dif==0)
    Index.use2 <- levelIndex2[use.ind2==1,]
    Index.use2$start <- c(1,cumsum(Index.use2$length)[-nrow(Index.use2)]+1)
    Index.use2$end <-  cumsum(Index.use2$length)
    
    coefs.noreg <- list()
    for(i in 1:nrow(Index.use2)){
        coefs.noreg[[i]] <- c(unlist(coefs.f[[i]]))
        names(coefs.noreg[[i]]) <- base.name[Index.use2$start[i]:Index.use2$end[i]]
    }
    
    output <- list("intercept"=intercept, "coefs"=coefs.noreg, "vcov"=vcov.f)
    
}



# ScreenINT ---------------------------------------------------------------

ScreenINT <- function(formula, data.main, data.x, screen, screen.type,
                      family, fac.level,screen.num.int, Gorder, maxIter,
                      internal.int, verbose=verbose){
    
    n.fac <- length(fac.level)
    ## verbose=FALSE, maxIter=maxIte
    term.f <- terms(formula, data=data.main)
    order.f <- attr(term.f, "order")
    all.vars0 <- c(all.vars(term.f)[1], attr(term.f, "term.labels")[order.f==1])
    if(screen==FALSE){        
        if(internal.int==FALSE){
            ## ############################ 
            ## Specified by Users.
            ## ############################ 
            ## Two-ways 
            if(any(order.f==2)==TRUE){
                indTwo.n <- strsplit(attr(term.f, "term.labels")[attr(term.f, "order")==2], "\\:")
                if(verbose==TRUE){
                    select.int <- as.data.frame(do.call("rbind", indTwo.n))
                    colnames(select.int) <- c("Fac.1","Fac.2")
                    cat("\nSpecified Two-way Interactions.\n")
                    print(select.int); rm(select.int)
                }
                indTwo <- matrix(NA, ncol=length(indTwo.n), nrow=2)
                for(j in 1:length(indTwo.n)){
                    indTwo[1,j] <- which(indTwo.n[[j]][1]==all.vars0) - 1
                    indTwo[2,j] <- which(indTwo.n[[j]][2]==all.vars0) - 1
                }
                ## indTwo <- apply(indTwo, 2, function(x) x[order(x)])
                ## indTwo <- indTwo[,order(indTwo[1,])]
                ## indTwo <- indTwo[,order(indTwo[2,])]
            }
            ## Three-ways 
            if(any(order.f==3)==TRUE){
                indThree.n <- strsplit(attr(term.f, "term.labels")[attr(term.f, "order")==3], "\\:")
                indThree <- matrix(NA, ncol=length(indThree.n), nrow=3)
                for(j in 1:length(indThree.n)){
                    indThree[1,j] <- which(indThree.n[[j]][1]==all.vars0) - 1
                    indThree[2,j] <- which(indThree.n[[j]][2]==all.vars0) - 1
                    indThree[3,j] <- which(indThree.n[[j]][3]==all.vars0) - 1
                }
                ## CHECK
                indThreeCheck <- indTwo2Three(indTwo=indTwo, n.fac=n.fac)
                if(is.null(indThreeCheck)){
                    check.three <- FALSE
                }else{
                    check.three <- c()
                    for(z in 1:ncol(indThree)){
                        check.three[z] <- any(apply(indThreeCheck == indThree[,j],2,all))
                    }
                }
                if(any(check.three==FALSE)){
                    stop("\nEvery pair of terms in Three-way Interactions should be specified in 'int2.formula' to satisfy strong hierarchy.\n")
                }
                rm(check.three); rm(indThreeCheck)
                
                if(verbose==TRUE){
                    select.int3 <- as.data.frame(do.call("rbind", indThree.n))
                    colnames(select.int3) <- c("Fac.1","Fac.2", "Fac.3")
                    cat("\nSpecified Three-way Interactions.\n")
                    print(select.int3); rm(select.int3)
                }
                Gorder <- 3
            }else{ indThree <- NULL}
        }else if(internal.int==TRUE & Gorder==2){
            ## ############################ 
            ## Analyze all Two-ways
            ## ############################ 
            if(verbose==TRUE) cat("\nAnalyzing all two-way interactions...\n")
            formula <- as.formula(paste(all.vars0[1]," ~ .*.", sep=""))
            indTwo <- combn(seq(1:n.fac),2)
            indThree <- NULL
        }else if(internal.int==TRUE & Gorder==3){
            ## ############################ 
            ## Analyze all Three-ways
            ## ############################ 
            if(verbose==TRUE) cat("\nAnalyzing all two-way and three-way interactions...\n")
            formula <- as.formula(paste(all.vars0[1]," ~ .*.*.", sep=""))
            indTwo <- combn(seq(1:n.fac),2)
            indThree <- combn(seq(1:n.fac),3)
        }
    }else if(screen==TRUE){
        ## ############################ 
        ## Data-driven Screening 
        ## ############################ 
        ## prepare the data for screening
        data.s <- matrix(NA, ncol=ncol(data.x),nrow=nrow(data.x))
        for(i in 1:ncol(data.s)){
            data.s[,i] <- as.numeric(data.x[,i]) - 1            
        }
        ## Do the screening here and modify formula.
        main.for <- paste(all.vars0[1],"~",paste(all.vars0[-1],collapse="+"),sep="")
        if(verbose==TRUE) cat("\n ***** Screening ***** \n")
        
        if(screen.type=="cv.min" | screen.type=="cv.1Std"){
            gl.screen <- glinternet.cv(X=data.s, Y=data.main[,1], family=family, numLevels=fac.level,
                                       nLambda=5, nFolds=10,verbose=FALSE, maxIter=maxIter)            
            if(screen.type=="cv.min"){
                if(is.null(gl.screen$activeSet[[1]]$catcat)==TRUE) indTwo <- NULL
                else indTwo <- t(gl.screen$activeSet[[1]]$catcat)
            }else if(screen.type=="cv.1Std"){
                if(is.null(gl.screen$activeSet1Std[[1]]$catcat)==TRUE) indTwo <- NULL
                else indTwo <- t(gl.screen$activeSet1Std[[1]]$catcat)
            }
            
        }else if (screen.type=="fixed"){
            gl.screen <- glinternet(X=data.s, Y=data.main[,1], family=family, numLevels=fac.level,
                                    numToFind=screen.num.int,
                                    verbose=FALSE, maxIter=maxIter)
            check.indTwo <- gl.screen$activeSet[[length(gl.screen$activeSet)]]$catcat
            if(is.null(check.indTwo)==FALSE){
                indTwo <- t(gl.screen$activeSet[[length(gl.screen$activeSet)]]$catcat)
            }else{
                indTwo <- NULL
            }
        }
        if(is.null(indTwo)==FALSE){
            if(verbose==TRUE){
                select.int <- as.data.frame(cbind(all.vars0[(indTwo[1,]+1)], all.vars0[(indTwo[2,]+1)]))
                colnames(select.int) <- c("Fac.1","Fac.2")
                cat("\nSelected Two-way Interactions.\n")
                print(select.int); rm(select.int)
            }
            
            int.for <- c()
            for(j in 1:ncol(indTwo)){
                int.for[j] <- paste(all.vars0[indTwo[1,j]+1], all.vars0[indTwo[2,j]+1], sep=":")
            }
            int.for <- paste(int.for, collapse="+")
        }else{
            if(verbose==TRUE) cat("\nNo Two-way Interaction is selected.\n")
            formula <- as.formula(main.for)
            indThree <- NULL
            Gorder <- 1
        }
        
        if(Gorder==3 & is.null(indTwo)==FALSE){
            indThree <- indTwo2Three(indTwo=indTwo, n.fac=n.fac)
            if(is.null(indThree)==FALSE){
                if(verbose==TRUE){
                    select.int3 <- as.data.frame(cbind(all.vars0[(indThree[1,]+1)],
                                                       all.vars0[(indThree[2,]+1)],
                                                       all.vars0[(indThree[3,]+1)]))
                    colnames(select.int3) <- c("Fac.1","Fac.2", "Fac.3")
                    cat("\nSelected Three-way Interactions.\n")
                    print(select.int3); rm(select.int3)
                }
                Gorder <- 3
                ## Formula
                int3.for <- c()
                for(j in 1:ncol(indThree)){
                    int3.for[j] <- paste(all.vars0[indThree[1,j]+1],
                                         all.vars0[indThree[2,j]+1],
                                         all.vars0[indThree[3,j]+1],sep=":")
                }
                int3.for <- paste(int3.for, collapse="+")
                formula <- as.formula(paste(main.for, int.for, int3.for, sep="+"))
            }else if(is.null(indTwo)==FALSE & is.null(indThree)==TRUE){
                if(verbose==TRUE) cat("\nNo Three-way Interaction is selected.\n")
                formula <- as.formula(paste(main.for, int.for, sep="+"))
                Gorder <- 2
            }
        }else if(Gorder==2 & is.null(indTwo)==FALSE){
            formula <- as.formula(paste(main.for, int.for, sep="+"))
            indThree <- NULL
            Gorder <- 2
        }
    }    
    output <- list("formula"=formula, "indTwo"=indTwo, "indThree"=indThree, "Gorder"=Gorder)
    return(output)
}


# cluster_se_glm ----------------------------------------------------------

cluster_se_glm <- function(model, cluster){
    
    #  Drop unused cluster indicators, if cluster var is a factor
    if (class(cluster) == "factor") {
        cluster <- droplevels(cluster)
    }
    
    if (nrow(model.matrix(model)) != length(cluster)) {
        stop("check your data: cluster variable has different N than model - you may have observations with missing data") 
    }
    
    M <- length(unique(cluster))
    N <- length(cluster)           
    K <- model$rank
    
    ## if(M<50) {
    ##     warning("Fewer than 50 clusters, variances may be unreliable (could try block bootstrap instead).")
    ## }
    
    dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
    uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
    rcse.cov <- dfc * sandwich(model, meat. = crossprod(uj)/N)
    return(rcse.cov)
}



# summary.CausaulANOVA ----------------------------------------------------

summary.CausalANOVA <-function(object, digit=4, verbose=TRUE, verbose.full=TRUE, ...){
    
    if(class(object)[2]=="stab"){
        fit <- object$fit
        stab.fit <- object$stab.fit
    }else{
        fit <- object
    }
    
    data <- fit$data
    Gorder <- fit$Gorder
    formula <- fit$formula
    formula.orig <- fit$formula.orig
    cost <- fit$cost
    n.fac <- length(fit$fac.level)
    fac.level <- fit$fac.level
    terms.f <- terms(formula, data=data)
    order.f <- attr(terms.f, "order")
    Name <- attr(terms.f, "term.labels")
    all.vars0 <- Name[order.f==1]
    inference <- fit$inference
    Stab <- as.numeric(class(object)[2]=="stab")
    level.list <- fit$level.list
    indTwo <- fit$indTwo
    indThree <- fit$indThree
    collapse <- fit$collapse
    screening <- fit$screen
    
    ## ################# 
    ## Compute Range
    ## ################# 
    range.reg <- c()
    for(z in 1:length(fit$coef)){
        range.reg[z] <-  max(fit$coef[[z]]) - min(fit$coef[[z]])
    }
    range <- round(range.reg,digits=3)
    
    Range <-as.data.frame(range)
    rownames(Range) <- Name
    
    ## Name Preparation 
    Factor.ame <- rep(Name[order.f==1], fac.level)
    Level.ame <- unlist(level.list)
    Base.ame <- unlist(lapply(fac.level, function(x) c(rep("",(x-1)), "***")))
    index.ame<- data.frame(cbind(Factor.ame, Level.ame, Base.ame))
    
    if(Gorder >=2){
        Name.amie2 <- Name[order.f==2]
        Factor.amie2 <- list()
        Level.amie2 <- list()
        Base.amie2 <- c()
        for(z in 1:sum(order.f==2)){
            size <- fac.level[[indTwo[1,z]]]*fac.level[[indTwo[2,z]]]
            Factor.amie2[[z]] <- rep(Name.amie2[z], times=size)
            Level.amie2[[z]] <- expand.grid(Var1=level.list[[indTwo[1,z]]], Var2=level.list[[indTwo[2,z]]])
            Base.amie2 <- c(Base.amie2, c(rep("", size-1), "***"))
        }
        Factor.amie2 <- unlist(Factor.amie2)
        Level.amie2 <- do.call("rbind",Level.amie2)
        index.amie2<- data.frame(cbind(Factor.amie2, Level.amie2, Base.amie2))
        
        
        if(Gorder == 3){
            Name.amie3 <- Name[order.f==3]
            Factor.amie3 <- list()
            Level.amie3 <- list()
            Base.amie3 <- c()
            for(z in 1:sum(order.f==3)){
                size <- fac.level[[indThree[1,z]]]*fac.level[[indThree[2,z]]]*fac.level[[indThree[3,z]]]
                Factor.amie3[[z]] <- rep(Name.amie3[z], times=size)
                Level.amie3[[z]] <- expand.grid(Var1=level.list[[indThree[1,z]]],
                                                Var2=level.list[[indThree[2,z]]],
                                                Var3=level.list[[indThree[3,z]]])
                Base.amie3 <- c(Base.amie3, c(rep("", size-1), "***"))
            }
            Factor.amie3 <- unlist(Factor.amie3)
            Level.amie3 <- do.call("rbind",Level.amie3)
            index.amie3<- data.frame(cbind(Factor.amie3, Level.amie3, Base.amie3))
        }                        
    }
    ## ################
    ## AME
    ## ################
    if(inference==TRUE){
        AME.CI <- fit$CI.table[order.f==1]
        AME.out <- do.call("rbind", AME.CI)
        AME.table <- round(AME.out, digit)
        AME.table2 <- data.frame(AME.table)
        
        AME.table.pr <- cbind(index.ame, AME.table2)
        colnames(AME.table.pr) <- c("Factor","Levels","base","AME","Std.Err","2.5%CI","97.5%CI")
        
        if(Gorder>=2){                        
            AMIE2.CI <- fit$CI.table[order.f==2]
            AMIE2.out <- do.call("rbind", AMIE2.CI)
            AMIE2.table <- round(AMIE2.out, digit)
            
            AMIE2.table.pr <- cbind(index.amie2, data.frame(AMIE2.table))
            colnames(AMIE2.table.pr) <- c("Factor","Level1","Level2", "base","AMIE","Std.Err","2.5%CI","97.5%CI")
            
            if(Gorder == 3){
                AMIE3.CI <- fit$CI.table[order.f==3]
                AMIE3.out <- do.call("rbind", AMIE3.CI)
                AMIE3.table <- round(AMIE3.out, digit)
                
                AMIE3.table.pr <- cbind(index.amie3, data.frame(AMIE3.table))
                colnames(AMIE3.table.pr) <- c("Factor","Level1","Level2","Level3",
                                              "base","AMIE","Std.Err","2.5%CI","97.5%CI")
            }else{
                AMIE3.table.pr <- NULL
            }
        }else{
            AMIE2.table.pr <- NULL
            AMIE3.table.pr <- NULL
        }
        
    }else if(inference==FALSE){
        AME.CI <- fit$coefs[order.f==1]
        AME.out0 <- list()
        for(z in 1:n.fac){
            AME.out0[[z]] <-  AME.CI[[z]] - AME.CI[[z]][length(AME.CI[[z]])]
        }
        AME.table <- data.frame(matrix(round(unlist(AME.out0),digit), ncol=1))
        
        AME.table.pr <- cbind(index.ame, AME.table)
        colnames(AME.table.pr) <- c("Factor","Levels","base","AME")
        
        if(Gorder>=2){
            AMIE2.CI <- fit$coefs[order.f==2]
            AMIE2.out0 <- list()
            for(z in 1:length(AMIE2.CI)){
                AMIE2.out0[[z]] <-  AMIE2.CI[[z]] - AMIE2.CI[[z]][length(AMIE2.CI[[z]])]
            }
            AMIE2.table <- data.frame(matrix(round(unlist(AMIE2.out0),digit), ncol=1))
            AMIE2.table.pr <- cbind(index.amie2, data.frame(AMIE2.table))
            colnames(AMIE2.table.pr) <- c("Factor","Level1","Level2", "base","AMIE")
            
            if(Gorder>=3){
                AMIE3.CI <- fit$coefs[order.f==3]
                AMIE3.out0 <- list()
                for(z in 1:length(AMIE3.CI)){
                    AMIE3.out0[[z]] <-  AMIE3.CI[[z]] - AMIE3.CI[[z]][length(AMIE3.CI[[z]])]
                }
                AMIE3.table <- data.frame(matrix(round(unlist(AMIE3.out0),digit), ncol=1))
                AMIE3.table.pr <- cbind(index.amie3, data.frame(AMIE3.table))
                colnames(AMIE3.table.pr) <- c("Factor","Level1","Level2","Level3", "base","AMIE")
            }else{
                AMIE3.table.pr <- NULL
            }           
        }else{
            AMIE2.table.pr <- NULL
            AMIE3.table.pr <- NULL
        }
        
        if(Stab==TRUE){
            ## Range 
            Range$select.prob <- round(stab.fit$range.stab[,2], digits=digit)
            colnames(Range)[2] <- "select.prob"
            
            ## AME
            AME.stab <- stab.fit$AME.stab[,2]      
            AME.table.pr$select.prob <- AME.stab            
        }
    }
    
    
    Range.pr <- as.data.frame(Range)
    
    ## Collapsing Results
    if(collapse==TRUE){
        collapse.u  <- Collapsing(fit)$collapse.level
        drop.fac <- unlist(lapply(collapse.u, function(x) all(x==1)))                        
        ## New Level Names
        collapse.name <- list()
        for(i in 1:n.fac){
            original.level <- level.list[[i]]
            new.name <- c()
            for(j in 1:length(unique(collapse.u[[i]]))){
                new.name[j] <- paste(original.level[collapse.u[[i]]==j],collapse="/")
            }
            collapse.name[[i]] <- unique(new.name[collapse.u[[i]]])
        }
        names(collapse.name) <- Name[order.f==1]
    }
    if(screening==TRUE){
        if(is.null(indTwo)==FALSE){
            select.int2 <- as.data.frame(cbind(all.vars0[(indTwo[1,])], all.vars0[(indTwo[2,])]))
            colnames(select.int2) <- c("Fac.1","Fac.2")
        }
        if(is.null(indTwo)==FALSE){
            select.int3 <- as.data.frame(cbind(all.vars0[(indThree[1,])], all.vars0[(indThree[2,])], all.vars0[(indThree[3,])]))
            colnames(select.int3) <- c("Fac.1","Fac.2", "Fac.3")
        }
    }
    
    
    if(verbose==TRUE){
        cat(" \nModel:\n ")
        print(formula)               
        
        cat("******\n")
        cat("Range:\n")
        cat("******\n")
        print(Range.pr)
        
        cat("*******************************\n")
        cat("Average Marginal Effects (AME):\n")
        cat("*******************************\n")
        print(AME.table.pr, row.names=F)
        cat("\n")
        
        if(verbose.full==TRUE & is.null(AMIE2.table.pr)==FALSE){
            cat("********************************************\n")
            cat("Average Marginal Interaction Effects (AMIE):\n")
            cat("********************************************\n")
            cat("Two-way Interactions:\n")
            print(AMIE2.table.pr, row.names=FALSE)
            
            if(is.null(AMIE3.table.pr)==FALSE){
                cat("\nThree-way Interactions:\n")
                print(AMIE3.table.pr, row.names=FALSE)
            }
        }
        
        if(screening==TRUE){
            cat("--------------------------------------------\n")
            cat("Selected Interactions; Results of Screening:\n")
            cat("--------------------------------------------\n")
            if(is.null(indTwo)==FALSE){
                print(select.int2)
            }
            if(is.null(indThree)==FALSE){
                print(select.int3)
            }
        }
        
        if(collapse==TRUE){
            cat("-----------------------------\n")
            cat("Results of Collapsing Levels:\n")
            cat("-----------------------------\n")
            print(collapse.name)           
        }
        
        if(inference==FALSE & Stab==FALSE){
            cat("===============================================================\n")
            cat(": To get confidence intervals, Use 'test.CausalANOVA' function:\n")
            cat("===============================================================\n")
        }
        
    }    
    output <- list("range"=Range.pr,"range.name"=Name,"AME"=AME.table.pr, "AMIE2"=AMIE2.table.pr,"AMIE3"=AMIE3.table.pr)
}




# test.CausaulANOVA -------------------------------------------------------


test.CausalANOVA <- function(fit, newdata, collapse.level=TRUE, diff=FALSE, pair.id=NULL,cluster=NULL){
    ## I need to install newdata as well as new information about pair.id
    
    eps <- fit$eps
    formula <- fit$formula
    main.formula <- fit$main.formula
    data <- model.frame(main.formula, fit$data)
    terms.f <- terms(formula, data=data)
    order.f <- attr(terms.f, "order")
    var.name <- attr(terms.f, "term.labels")[order.f==1]
    fac.level <- fit$fac.level
    n.fac <- length(fac.level)
    data <- fit$data
    family <- fit$family
    
    data.main0 <- model.frame(main.formula, data=newdata)
    
    if(collapse.level==TRUE){
        collapse.u <- collapse <- Collapsing(fit)$collapse.level
        drop.fac <- unlist(lapply(collapse, function(x) all(x==1)))
        
        ## Recover Collapsed Factor
        for(z in which(drop.fac)){
            n1 <- ceiling(length(collapse.u[[z]])/2)
            n2 <- length(collapse.u[[z]]) - n1
            collapse.u[[z]] <- c(rep(1, n1), rep(2,n2))
        }
        
        ## New Level Names
        collapse.u2 <- list()
        for(i in 1:n.fac){
            original.level <- levels(data[,(i+1)])
            new.name <- c()
            for(j in 1:length(unique(collapse.u[[i]]))){
                new.name[j] <- paste(original.level[collapse.u[[i]]==j],collapse="/")
            }
            collapse.u2[[i]] <- new.name[collapse.u[[i]]]
        }
        
        ## ############################ 
        ## Reconstruct Formula
        ## ############################ 
        
        fac.table <- attr(terms.f, "factors")[-1,]
        d.table <- fac.table[drop.fac==TRUE,]
        
        if(sum(drop.fac)==0) drop.ind <- 0
        if(sum(drop.fac)==1) drop.ind <- (d.table > 0)*(order.f>1)
        if(sum(drop.fac)>1){
            drop.ind <- (apply(d.table,2,sum) > 0)*(order.f>1)
        }
        
        term.labels <- attr(terms.f, "term.labels")
        
        if(all(drop.fac==FALSE) | all(drop.ind==0)==TRUE){
            new.formula <- formula
        }else{
            new.formula <-
                as.formula(paste(as.character(formula)[2], "~", paste(term.labels[drop.ind==0], collapse="+"),sep=""))
        }
    }else{
        new.formula <- formula
    }
    ## ############################
    ## Prepare New Data 
    ## ############################
    data.main <- model.frame(new.formula, data=data.main0)
    data.x <- model.frame(new.formula, data=data.main0)[,-1]
    terms.f.new <- terms(new.formula, data=data.main0)
    nway.new <- max(attr(terms.f.new, "order"))
    ord.fac <- rep(FALSE, n.fac)
    
    ## Releveling New Data
    if(collapse.level==TRUE){
        for(i in 1:n.fac){
            levels(data.main[,(i+1)]) <- collapse.u2[[i]]
            data.main[,(i+1)] <- factor(data.main[,(i+1)],ordered=FALSE)
        }
    }
    
    fit.new <- CausalANOVAFit(formula = new.formula, internal.int=FALSE,
                              pair.id=pair.id, diff=diff, 
                              data=data.main, cost=1, cluster=cluster,
                              family=family, nway=nway.new,
                              eps=eps, screen=FALSE, ord.fac=NULL,
                              collapse=FALSE, verbose=FALSE)
    
    coefs <- fit.new$coefs
    vcov  <- fit.new$vcov
    fit.new$inference <- TRUE
    
    ## Index
    levelIndex <- CreatelevelIndex(fac.level=fit.new$fac.level, ord.fac=fit.new$ord.fac, Gorder=fit.new$Gorder,
                                   indTwo=fit.new$indTwo, indThree=fit.new$indThree)
    use.ind <-  (levelIndex$plus==1)*(levelIndex$dif==0)
    Index.use <- levelIndex[use.ind==1,]
    Index.use$start <- c(1,cumsum(Index.use$length)[-nrow(Index.use)]+1)
    Index.use$end <-  cumsum(Index.use$length)
    
    order.f <- attr(terms.f.new, "order")
    
    result.tab <- list()
    for(z in 1:length(coefs)){
        seq.u <- Index.use$start[z]:Index.use$end[z]
        effect <- coefs[[z]] - coefs[[z]][length(coefs[[z]])]
        vcov.p <- c()
        for(i in 1:(length(seq.u)-1)){
            vcov.p[i] <- VarEffect(ind1=seq.u[i],ind2=seq.u[length(seq.u)],vcov.full=vcov)
        }
        sd <- c(sqrt(vcov.p),0)
        table <- as.data.frame(cbind(effect, sd, effect - 1.96*sd, effect + 1.96*sd))
        if(Index.use$order[z]==1){
            colnames(table) <- c("AME", "sd", "2.5% CI", "97.5% CI")
        }else{
            colnames(table) <- c("AMIE", "sd", "2.5% CI", "97.5% CI")
        }
        result.tab[[z]] <- table
    }
    
    fit.new[[length(fit.new)+1]] <- result.tab
    names(fit.new)[length(fit.new)] <- "CI.table"
    ## output <- c(fit.new, "result.tab"=result.tab)
    
    output <- fit.new
    class(output) <- c("CausalANOVA","list")
    return(output)
}


# Lcombinefunction --------------------------------------------------------

Lcombinefunction <- function(fac.level,ord.fac,facCons,Gorder=3, indTwo=NULL, indThree=NULL){
    
    ## Make the difference in three-way interactions
    D3function <- function(alevel,blevel,clevel,type="A",ordABC=c(FALSE,FALSE,FALSE)){
        if(type=="A"){
            D3 <- D3Afunction(alevel=alevel,blevel=blevel,clevel=clevel,ord=ordABC[1])
        }else if(type=="B"){
            D3 <- D3Bfunction(alevel=alevel,blevel=blevel,clevel=clevel,ord=ordABC[2])
        }else if(type=="C"){
            D3 <- D3Cfunction(alevel=alevel,blevel=blevel,clevel=clevel,ord=ordABC[3])
        }
        return(D3)
    }
    
    D3Afunction <- function(alevel,blevel,clevel,ord=FALSE){
        D3a <- matrix(0,nrow=0,ncol=(alevel*blevel*clevel))  
        if(ord==FALSE){
            for(i in 1:(alevel-1)){
                for(j  in (i+1):alevel){
                    base.v <- rep(0,alevel)
                    base.v[i] <- -1
                    base.v[j] <- 1
                    base.mat <- diag((blevel*clevel)) %x% t(base.v)
                    D3a <- rbind(D3a,base.mat)
                }
            }
        }else if(ord==TRUE){
            for(i in 1:(alevel-1)){
                base.v <- rep(0,alevel)
                base.v[i] <- -1
                base.v[(i+1)] <- 1
                base.mat <- diag((blevel*clevel)) %x% t(base.v)
                D3a <- rbind(D3a,base.mat)
            }        
        }
        return(D3a)
    }
    
    ## (B): Make the difference in two-way interactions 
    D3Bfunction <- function(alevel,blevel,clevel,ord=FALSE){    
        if(ord==FALSE){
            D3b <- matrix(0,nrow=0,ncol=(alevel*blevel*clevel))
            for(i in 1:(blevel-1)){     
                left.temp <- matrix(0,ncol=((i-1)*(alevel)),nrow=(alevel))
                left.main <- -diag((alevel))
                left.final <- cbind(left.temp,left.main)
                for(j  in (i+1):blevel){
                    middle.temp <-  matrix(0,ncol=((j-i-1)*(alevel)),nrow=(alevel))
                    right.main <- diag((alevel))
                    right.temp <- matrix(0,ncol=((blevel-j)*(alevel)),nrow=(alevel))
                    right.final <- cbind(middle.temp,right.main,right.temp)
                    D2b <- cbind(left.final,right.final) ## The base
                    D3base <- diag(clevel) %x% D2b
                    D3b <- rbind(D3b, D3base)
                }
            }
        }else if (ord==TRUE){
            D3b <- matrix(0,nrow=0,ncol=(alevel*blevel*clevel))
            for(i in 1:(blevel-1)){     
                left <- matrix(0,ncol=((i-1)*(alevel)),nrow=(alevel))
                Main <- cbind(-diag((alevel)),diag((alevel)))
                right <- matrix(0,ncol=((blevel-(i+1))*(alevel)),nrow=(alevel))
                D2b <- cbind(left,Main,right)
                D3base <- diag(clevel) %x% D2b
                D3b <- rbind(D3b, D3base)
            }
        }
        return(D3b)
        ## psy_B(2,3), psy_B(2,4), psy_B(2,5) 
    }
    
    D3Cfunction <- function(alevel,blevel,clevel,ord=FALSE){
        if(ord==FALSE){
            D3c <- matrix(0,nrow=0,ncol=(alevel*blevel*clevel))
            for(i in 1:(clevel-1)){     
                left.temp <- matrix(0,ncol=((i-1)*(alevel*blevel)),nrow=(alevel*blevel))
                left.main <- -diag((alevel*blevel))
                left.final <- cbind(left.temp,left.main)
                for(j  in (i+1):clevel){
                    middle.temp <-  matrix(0,ncol=((j-i-1)*(alevel*blevel)),nrow=(alevel*blevel))
                    right.main <- diag((alevel*blevel))
                    right.temp <- matrix(0,ncol=((clevel-j)*(alevel*blevel)),
                                         nrow=(alevel*blevel))
                    right.final <- cbind(middle.temp,right.main,right.temp)
                    D3c <- rbind(D3c,cbind(left.final,right.final))
                }
            }
        }else if (ord==TRUE){
            Left <- cbind(-diag((alevel*blevel)*(clevel-1)),
                          matrix(0,ncol=(alevel*blevel),nrow=(alevel*blevel)*(clevel-1)))
            Right <- cbind(matrix(0,ncol=(alevel*blevel),nrow=(alevel*blevel)*(clevel-1)),
                           diag((alevel*blevel)*(clevel-1)))
            D3c <- Left + Right
        }
        return(D3c)
    }
    
    
    ## From this version, I use D2 rather than D3 
    ## Make the difference in two-way interactions 
    D2function <- function(alevel,blevel,type,ordAB){
        if(type=="A"){
            D2 <- D2Afunction(alevel=alevel,blevel=blevel,ord=ordAB[1])
        }else if(type=="B"){        
            D2 <- D2Bfunction(alevel=alevel,blevel=blevel,ord=ordAB[2])
        }
        return(D2)
        ## D2 matrix
        ## nrow = differences between two-way interactions fixing the sublevel
        ## ncol = alevel*blevel (coefficients for two-ways)
    }
    
    ## (A): Make the difference in two-way interactions 
    D2Afunction <- function(alevel,blevel,ord=FALSE){
        D2a <- matrix(0,nrow=0,ncol=(alevel*blevel))
        if(ord==FALSE){
            for(i in 1:(alevel-1)){
                for(j  in (i+1):alevel){
                    base.v <- rep(0,alevel)
                    base.v[i] <- -1
                    base.v[j] <- 1
                    base.mat <- diag(blevel) %x% t(base.v)
                    D2a <- rbind(D2a,base.mat)
                }
            }
        }else if(ord==TRUE){
            for(i in 1:(alevel-1)){
                base.v <- rep(0,alevel)
                base.v[i] <- -1
                base.v[(i+1)] <- 1
                base.mat <- diag(blevel) %x% t(base.v)
                D2a <- rbind(D2a,base.mat)
            }        
        }
        return(D2a)
    }
    
    ## (B): Make the difference in two-way interactions 
    D2Bfunction <- function(alevel,blevel,ord=FALSE){    
        if(ord==FALSE){
            D2b <- matrix(0,nrow=0,ncol=(alevel*blevel))
            for(i in 1:(blevel-1)){     
                left.temp <- matrix(0,ncol=((i-1)*(alevel)),nrow=(alevel))
                left.main <- -diag((alevel))
                left.final <- cbind(left.temp,left.main)
                for(j  in (i+1):blevel){
                    middle.temp <-  matrix(0,ncol=((j-i-1)*(alevel)),nrow=(alevel))
                    right.main <- diag((alevel))
                    right.temp <- matrix(0,ncol=((blevel-j)*(alevel)),nrow=(alevel))
                    right.final <- cbind(middle.temp,right.main,right.temp)
                    D2b <- rbind(D2b,cbind(left.final,right.final))
                }
            }
        }else if (ord==TRUE){
            Left <- cbind(-diag((alevel)*(blevel-1)),matrix(0,ncol=alevel,nrow=(alevel)*(blevel-1)))
            Right <- cbind(matrix(0,ncol=alevel,nrow=(alevel)*(blevel-1)), diag((alevel)*(blevel-1)))
            D2b <- Left + Right
        }
        return(D2b)
        ## psy_B(2,3), psy_B(2,4), psy_B(2,5) 
    }
    
    ## (Main): Making the difference between all levels
    D1function <- function(nlevel,ord=FALSE){
        if(ord==FALSE){
            qp = (nlevel)*(nlevel-1)/2
            D1.mat <- matrix(0,ncol=(nlevel),nrow=0)    
            if (qp==1) D1.mat <- matrix(c(-1,1),ncol=2,nrow=1)
            if (qp>1){
                for(i in 1 : (nlevel-1)){
                    w.diag <- diag((nlevel-i))
                    left.v <- c(rep(0,(i-1)),-1)
                    left.mat <- matrix(rep(left.v,(nlevel-i)),nrow=(nlevel-i),byrow=TRUE)
                    w.mat <- cbind(left.mat,w.diag)
                    D1.mat <- rbind(D1.mat,w.mat)
                }
            }
        }else if(ord==TRUE){
            if(nlevel==2){
                D1.mat <- matrix(c(-1,1),ncol=2,nrow=1)
            }else{
                D1.mat <- cbind(0,diag(nlevel-1)) + cbind(-diag(nlevel-1),0)
            }
            
        }    
        return(D1.mat)
        ## nrow = differences between main effects
        ## ncol = mainlevel (coefficients for the main factor)
    }
    
    L1function <- function(nlevel,ord=FALSE){
        ## Creating L for the main effects
        D1 <- D1function(nlevel,ord=ord)
        if(ord==FALSE){
            qp <- (nlevel)*(nlevel-1)/2
        }else if(ord==TRUE){
            qp <- nlevel-1
        }
        L1 <- cbind(D1, -D1, -diag(qp), diag(qp))
        return(L1)
    }
    
    ## From this version, we call this L2 function.
    
    L2function <- function(alevel,blevel,ordAB=c(FALSE,FALSE)){
        ## Creating L for interaction effects
        if(ordAB[1]==FALSE){
            r1 <- (blevel)*(alevel)*(alevel-1)/2
        }else if(ordAB[1]==TRUE){
            r1 <- (blevel)*(alevel-1)
        }
        if(ordAB[2]==FALSE){
            r2 <- (alevel)*(blevel)*(blevel-1)/2
        }else if(ordAB[2]==TRUE){
            r2 <- (alevel)*(blevel-1)
        }
        D2a <- D2function(alevel=alevel,blevel=blevel,type="A",ordAB=ordAB)
        D2b <- D2function(alevel=alevel,blevel=blevel,type="B",ordAB=ordAB)    
        zero.r1 <- matrix(0,ncol=r1,nrow=r2)
        zero.r2 <- matrix(0,ncol=r2,nrow=r1)
        L2top <- cbind(D2a, -D2a, -diag(r1), diag(r1), zero.r2, zero.r2)
        L2bottom <- cbind(D2b, -D2b, zero.r1, zero.r1, -diag(r2), diag(r2))
        L2 <- rbind(L2top,L2bottom)
        return(L2)
    }
    
    L3function <- function(alevel,blevel,clevel,ordABC=c(FALSE,FALSE,FALSE)){
        ## Creating L for interaction effects
        if(ordABC[1]==FALSE){
            r1 <- (blevel*clevel)*(alevel)*(alevel-1)/2
        }else if(ordABC[1]==TRUE){
            r1 <- (blevel*clevel)*(alevel-1)
        }
        if(ordABC[2]==FALSE){
            r2 <- (alevel*clevel)*(blevel)*(blevel-1)/2
        }else if(ordABC[2]==TRUE){
            r2 <- (alevel*clevel)*(blevel-1)
        }
        if(ordABC[3]==FALSE){
            r3 <- (alevel*blevel)*(clevel)*(clevel-1)/2
        }else if(ordABC[3]==TRUE){
            r3 <- (alevel*blevel)*(clevel-1)
        }
        D3a <- D3function(alevel=alevel,blevel=blevel,clevel=clevel,type="A",ordABC=ordABC)
        D3b <- D3function(alevel=alevel,blevel=blevel,clevel=clevel,type="B",ordABC=ordABC)
        D3c <- D3function(alevel=alevel,blevel=blevel,clevel=clevel,type="C",ordABC=ordABC)
        zero.1 <- matrix(0,ncol=(2*(r2+r3)),nrow=r1)
        zero.2l <- matrix(0,ncol=(2*r1),nrow=r2)
        zero.2r <- matrix(0,ncol=(2*r3),nrow=r2)
        zero.3 <- matrix(0,ncol=(2*(r1+r2)),nrow=r3)
        L3.1 <- cbind(D3a, -D3a, -diag(r1), diag(r1), zero.1)
        L3.2 <- cbind(D3b, -D3b, zero.2l, -diag(r2), diag(r2), zero.2r)
        L3.3 <- cbind(D3c, -D3c, zero.3, -diag(r3), diag(r3))
        L3 <- rbind(L3.1,L3.2,L3.3)
        return(L3)
    }
    
    
    
    ## fac.level <- c(4,7,4)
    ## ord <- c(TRUE,TRUE,TRUE)
    
    ## maybe we can change here. 
    levelIndex <- CreatelevelIndex(fac.level=fac.level,ord.fac=ord.fac,Gorder=Gorder,
                                   indTwo=indTwo, indThree=indThree)
    use.ind <- (levelIndex$plus==1) * (levelIndex$dif==0)
    Index.use <- levelIndex[use.ind==1,]
    Fac.index <- Index.use[,regexpr("Fac",colnames(Index.use))>0]
    coef.length <- sum(levelIndex$length)
    
    
    L.list <- list()
    
    for(z in 1:nrow(Index.use)){
        if(Index.use$order[z]==1){
            level.main <- fac.level[which(Fac.index[z,]==1)]
            ord.main <- ord.fac[which(Fac.index[z,]==1)]
            L.list[[z]] <- L1function(nlevel=level.main,ord=ord.main)
        }else if(Index.use$order[z]==2){
            level.main <- fac.level[which(Fac.index[z,]==1)]
            ord.main <- ord.fac[which(Fac.index[z,]==1)]
            L.list[[z]] <- L2function(alevel=level.main[1],
                                      blevel=level.main[2],ordAB=ord.main)
        }else if(Index.use$order[z]==3){
            level.main <- fac.level[which(Fac.index[z,]==1)]
            ord.main <- ord.fac[which(Fac.index[z,]==1)]
            L.list[[z]] <- L3function(alevel=level.main[1],blevel=level.main[2],
                                      clevel=level.main[3],
                                      ordABC=ord.main)
        }
    }   
    
    Lm <- matrix(0,nrow=0,ncol=(coef.length))
    ## Need to make this generalizable
    for(l in 1:length(L.list)){
        L.minilist <- list()        
        for(ind in 1:length(L.list)){
            if(ind!=l){
                L.minilist[[ind]] <- matrix(0,nrow=nrow(L.list[[l]]),
                                            ncol=ncol(L.list[[ind]]))
            }else if(ind==l){
                L.minilist[[ind]] <- L.list[[l]]
            }            
        }
        Lb <- do.call(cbind,L.minilist)
        Lm <- rbind(Lm,Lb)
    }
    
    l.slack <- lengthSlack(fac.level=fac.level,ord.fac=ord.fac,Gorder=Gorder,
                           facCons=facCons)$l.slack
    Lm <- cbind(0,Lm)
    Lm <- rbind(0, Lm)
    Lslack <- matrix(0,nrow=nrow(Lm),ncol=l.slack)
    
    L <- cbind(Lm,Lslack)
    
    return(L)
}

# lengthSlack -------------------------------------------------------------

lengthSlack <- function(fac.level,ord.fac,Gorder=2, facCons=FALSE){
    ## This does not depend on Gorder.
    if(facCons==FALSE){
        n.fac <- length(fac.level)
        length.slack <- rep(0,times=n.fac)
        
        for(i in 1:n.fac){
            if(ord.fac[i]==FALSE){
                length.slack[i] <- fac.level[i]*(fac.level[i]-1)/2
            }else{
                length.slack[i] <- (fac.level[i]-1)
            }
        }
        l.slack <- sum(length.slack)
        output <- list("l.slack"=l.slack,
                       "length.full"=length.slack)
    }else if(facCons==TRUE){
        n.fac <- length(fac.level)
        first <- n.fac
        second <- choose(n.fac,2)
        if(Gorder==3){
            third <- choose(n.fac,3)
            l.slack <- first + second + third
        }else if(Gorder==2){
            l.slack <- first + second 
        }
        output <- list("l.slack"=l.slack)
    }
    return(output)
}

# PsyConstraintCombine ----------------------------------------------------

PsyConstraintCombine <- function(fac.level,ord.fac, Gorder, indTwo=NULL, indThree=NULL){
    
    PsyIJfunction <- function(i,j,type.ind,plus,fac.level,ord.fac,Gorder, indTwo, indThree){
        plus.t <- ifelse(plus==TRUE,1,0)
        
        ## Actually, level.vec will do the most work to keep track of names
        levelIndex <- CreatelevelIndex(fac.level=fac.level,ord.fac=ord.fac,
                                       Gorder=Gorder, indTwo=indTwo, indThree=indThree)
        Fac.index <- levelIndex[,regexpr("Fac",colnames(levelIndex))>0]
        use.ind <-  (levelIndex$plus==plus.t)*(levelIndex$dif==1)*(Fac.index[,type.ind]==2)
        ## Here we keep track of all factors and factor interactions related to "type.ind". 
        Index.use <- levelIndex[use.ind==1,]
        Fac.Ind.use <- Fac.index[use.ind==1,]
        Index.ind <- seq(1:nrow(levelIndex))[use.ind==1]
        
        mainlevel <- fac.level[type.ind]
        ordMain <- ord.fac[type.ind]
        
        ## Create the Basline Vectors
        if(ordMain==FALSE){
            level.comb <- combn(seq(from=1,to=mainlevel),2)
            ind.c <- which(((level.comb[1,]==i) * (level.comb[2,]==j))==1)
            base <- rep(0,times=((mainlevel)*(mainlevel-1)/2))
            base[ind.c] <- 1
        }else if(ordMain==TRUE){
            base <- rep(0,times=(mainlevel-1))
            base[i] <- 1
        }
        
        ## First Row (Main Effects)
        left.n.col1 <- sum(levelIndex$length[1:(Index.ind[1]-1)])
        right.n.col1 <- sum(levelIndex$length[(Index.ind[1]+1):nrow(levelIndex)])
        if(is.na(right.n.col1)==FALSE){
            PsyIJ <- c(rep(0,left.n.col1),base,rep(0,right.n.col1))
        }else{
            PsyIJ <- c(rep(0,left.n.col1),base)
        }
        
        ## ############################ 
        ## For Two-way and Three-way
        ## ############################ 
        ## Create Sublevel Matrix
        if(nrow(Index.use)>=2){
            Fac.sub <-t(t(Fac.Ind.use==1) * fac.level)
            sublevel <- rep(0,times=nrow(Fac.Ind.use))
            
            for(z in 2:nrow(Index.use)){
                sublevel[z] <- prod(Fac.sub[z,][Fac.sub[z,]!=0])
            }
            
            ## The base matrix for each set of coefficients
            Base.P.list <- list()
            for(z in 1:nrow(Index.use)){
                if(z>1){
                    Base.P.list[[z]] <- t(base) %x% diag(sublevel[z])
                }else{
                    Base.P.list[[z]] <- base
                }
            }
            
            
            ## Gather (Two way Effects)
            ## Need to make this generalizable
            ## I need to refer to the LevelIndex 
            for(p in 2:length(Base.P.list)){
                left.n.col <- sum(levelIndex$length[1:(Index.ind[p]-1)])
                Pb.left <- matrix(0,ncol=left.n.col,nrow=nrow(Base.P.list[[p]]))
                if((nrow(levelIndex)-Index.ind[p])>=2){
                    right.n.col <- sum(levelIndex$length[(Index.ind[p]+1):nrow(levelIndex)])         
                    Pb.right <- matrix(0,ncol=right.n.col,nrow=nrow(Base.P.list[[p]]))
                    Pb <- cbind(Pb.left, Base.P.list[[p]],Pb.right)
                }else if(nrow(levelIndex)==(1+Index.ind[p])){
                    Pb <- cbind(Pb.left, Base.P.list[[p]],
                                matrix(0,ncol=ncol(Base.P.list[[p]]),nrow=nrow(Base.P.list[[p]])))
                }else if(nrow(levelIndex)==(Index.ind[p])){
                    Pb <- cbind(Pb.left, Base.P.list[[p]])
                }
                PsyIJ <- rbind(PsyIJ,Pb)
            }
            ## Add a column for Mu.
            PsyIJ <- cbind(0, PsyIJ)
        }else{
            ## No interaction case
            ## Add a column for Mu.
            PsyIJ <- matrix(c(0, PsyIJ), nrow=1)
        }
        return(PsyIJ)
    }
    
    ## a <- PsyConstraint(type.ind=3,fac.level,ord.fac)
    
    
    
    ## Turn PhyMatrix Intro Constraint Form
    PsyConstraintFunction <- function(type.ind,fac.level,ord.fac,Gorder,
                                      indTwo, indThree){
        ## Main Job is to combine PsyIJ and Matrix for slack variables
        
        levelIndex <- CreatelevelIndex(fac.level,ord.fac, Gorder,indTwo=indTwo, indThree=indThree)
        coef.length <- sum(levelIndex$length)
        coef.length.mu <- coef.length + 1
        
        l.slack <- lengthSlack(fac.level=fac.level,ord.fac=ord.fac)$l.slack
        slack.full <- lengthSlack(fac.level=fac.level,ord.fac=ord.fac)$length.full
        
        mainlevel <- fac.level[type.ind]
        ordMain <- ord.fac[type.ind]
        if(type.ind>1){
            jump <- cumsum(slack.full)[(type.ind-1)]
        }else{
            jump <- 0
        }
        
        PsyFinal <- matrix(0,nrow=0,ncol=(coef.length.mu+l.slack))
        
        if(ordMain==FALSE){
            for(i in 1:(mainlevel-1)){    
                for(j in (i+1):mainlevel){
                    PijPlus <- PsyIJfunction(i=i,j=j,type.ind=type.ind,plus=TRUE,
                                             fac.level=fac.level,ord.fac=ord.fac,
                                             Gorder=Gorder, indTwo=indTwo, indThree=indThree)
                    PijMinus <- PsyIJfunction(i=i,j=j,type.ind=type.ind,plus=FALSE,
                                              fac.level=fac.level,ord.fac=ord.fac,
                                              Gorder=Gorder,indTwo=indTwo, indThree=indThree)                
                    Pij <- PijPlus + PijMinus
                    
                    ## Create Matrix for Slack variables
                    base <- rep(0,times=l.slack)
                    level.comb <- combn(seq(from=1,to=mainlevel),2)
                    ind <- which(((level.comb[1,]==i) * (level.comb[2,]==j))==1)
                    ind.use <- ind + jump
                    base[ind.use] <- 1
                    PsySlackIJ <- t(base %x% t(rep(1,times=nrow(Pij))))
                    PsyFinalIJ <- cbind(Pij, -PsySlackIJ)
                    PsyFinal <- rbind(PsyFinal,PsyFinalIJ)
                }
            }
        }else if(ordMain==TRUE){
            for(i in 1:(mainlevel-1)){    
                PijPlus <- PsyIJfunction(i=i,j=(i+1),type.ind=type.ind,plus=TRUE,
                                         fac.level=fac.level,ord.fac=ord.fac,
                                         Gorder=Gorder,indTwo=indTwo, indThree=indThree)
                PijMinus <- PsyIJfunction(i=i,j=(i+1),type.ind=type.ind,plus=FALSE,
                                          fac.level=fac.level,ord.fac=ord.fac,
                                          Gorder=Gorder,indTwo=indTwo, indThree=indThree)
                Pij <- PijPlus + PijMinus
                
                ## Create Matrix for Slack variables
                base <- rep(0,times=l.slack)
                ind <- i
                ind.use <- ind + jump
                base[ind.use] <- 1
                PsySlackIJ <- t(base %x% t(rep(1,times=nrow(Pij))))                      
                PsyFinalIJ <- cbind(Pij, -PsySlackIJ)
                PsyFinal <- rbind(PsyFinal,PsyFinalIJ)
            }        
        }    
        return(PsyFinal)
    }
    
    
    n.fac <- length(fac.level)
    levelIndex <- CreatelevelIndex(fac.level,ord.fac,Gorder,indTwo=indTwo, indThree=indThree)
    coef.length <- sum(levelIndex$length)
    coef.length.mu <- coef.length + 1    
    l.slack <- lengthSlack(fac.level=fac.level,ord.fac=ord.fac)$l.slack
    PsyConstraintF <- matrix(0, ncol=(coef.length.mu+l.slack),nrow=0)
    for(psy in 1:n.fac){
        Psy.temp <- PsyConstraintFunction(type.ind=psy,fac.level=fac.level,ord.fac=ord.fac,
                                          Gorder=Gorder, indTwo=indTwo,indThree=indThree)
        PsyConstraintF <- rbind(PsyConstraintF, Psy.temp)
    }
    return(PsyConstraintF)
}

# CreateANOVAconst -------------------------------------------------------
CreateANOVAconst <- function(fac.level,ord.fac,Gorder,facCons, indTwo=NULL, indThree=NULL){
    ## Two-way Interactions
    TwoWayConst <- function(alevel,blevel,type="alpha"){
        ## type determines which to fix
        if(type=="beta"){
            ## Fix j (beta for AB)
            int.constAB <- diag(blevel) %x% t(rep(1,times=alevel))
        }else if(type=="alpha"){
            ## Fix i (alpha for AB)
            int.constAB <- do.call(cbind, replicate(blevel, diag(alevel), simplify=FALSE))
        }
        return(int.constAB)
    }
    
    ThreeWayConst <- function(alevel,blevel,clevel,type="alpha"){
        ## type determines which to average over
        if(type=="gamma"){
            ## Fix j (beta for AB)
            int.constABC <- do.call(cbind, replicate(clevel, diag(alevel*blevel), 
                                                     simplify=FALSE))
        }else if(type=="beta"){
            int.constBase <- do.call(cbind, replicate(blevel, diag(alevel), 
                                                      simplify=FALSE))
            ## Fix i (alpha for AB)
            int.constABC <- diag(clevel) %x% int.constBase
        }else if(type=="alpha"){
            int.constABC <- diag(blevel*clevel) %x% t(rep(1,times=(alevel)))
        }
        return(int.constABC)
    }
    
    
    ## I introduce Slack Matrix later
    ## I will add intercept later
    
    ## Even with ordered Categorical varialbe, almost the same
    
    ## We can fix this. 
    levelIndex <- CreatelevelIndex(fac.level=fac.level,ord.fac=ord.fac,Gorder=Gorder,
                                   indTwo=indTwo, indThree=indThree)
    Fac.index <- levelIndex[,regexpr("Fac",colnames(levelIndex))>0]
    use.ind <-  (levelIndex$plus==1)*(levelIndex$dif==0)
    Index.use <- levelIndex[use.ind==1,]
    Fac.Ind.use <- Fac.index[use.ind==1,]
    Index.ind <- seq(1:nrow(levelIndex))[use.ind==1]
    coef.length <- sum(levelIndex$length)
    l.slack <- lengthSlack(fac.level=fac.level,ord.fac=ord.fac,Gorder=Gorder,
                           facCons=facCons)$l.slack       
    
    ## The base matrix for each set of coefficients
    Base.C.list <- list()
    for(z in 1:nrow(Index.use)){
        if(Index.use$order[z]==1){
            mainlevel <- fac.level[z]
            Base.C.list[[z]] <- rep(1,mainlevel)
        }else if(Index.use$order[z]==2){
            level.temp <- (Fac.Ind.use[z,]==1) * fac.level
            mainlevel <- level.temp[level.temp!=0]
            ## for the same interaction coefficients, so we can take rbind here. 
            Base.C.list[[z]] <- rbind(TwoWayConst(alevel=mainlevel[1],blevel=mainlevel[2],type="alpha"),
                                      TwoWayConst(alevel=mainlevel[1],blevel=mainlevel[2],type="beta"))
            
        }else if(Index.use$order[z]==3){
            level.temp <- (Fac.Ind.use[z,]==1) * fac.level
            mainlevel <- level.temp[level.temp!=0]           
            Base.C.list[[z]] <- rbind(ThreeWayConst(alevel=mainlevel[1],blevel=mainlevel[2],
                                                    clevel=mainlevel[3],type="alpha"),
                                      ThreeWayConst(alevel=mainlevel[1],blevel=mainlevel[2],
                                                    clevel=mainlevel[3],type="beta"),
                                      ThreeWayConst(alevel=mainlevel[1],blevel=mainlevel[2],
                                                    clevel=mainlevel[3],type="gamma"))
        }
    }
    
    ## Combine into the larger Matrix using the index 
    Const <- matrix(0,nrow=0,ncol=coef.length)
    for(z in 1:nrow(Index.use)){
        if(Index.use$order[z]==1){
            right.n.col1 <- sum(levelIndex$length[(Index.ind[z]+2):nrow(levelIndex)])
            if(Index.ind[z]>1){
                left.n.col1 <- sum(levelIndex$length[1:(Index.ind[z]-1)])
                Const.temp <- c(rep(0,left.n.col1),Base.C.list[[z]],-Base.C.list[[z]], rep(0,right.n.col1))
            }else{
                Const.temp <- c(Base.C.list[[z]],-Base.C.list[[z]], rep(0,right.n.col1))
            }
            Const <- rbind(Const, Const.temp)
        }else if(Index.use$order[z]>1){
            left.n.col <- sum(levelIndex$length[1:(Index.ind[z]-1)])
            right.n.col <- sum(levelIndex$length[(Index.ind[z]+2):nrow(levelIndex)])
            Cb.left <- matrix(0,ncol=left.n.col,nrow=nrow(Base.C.list[[z]]))
            Cb.right <- matrix(0,ncol=right.n.col,nrow=nrow(Base.C.list[[z]]))
            Const.temp <- cbind(Cb.left, Base.C.list[[z]], -Base.C.list[[z]], Cb.right)
            Const <- rbind(Const, Const.temp)
        }
    }
    
    ## Add intercept and slack variables
    ConstFinal <- cbind(0, Const, matrix(0,nrow=nrow(Const),ncol=l.slack))
    return(ConstFinal)
}


# Zcombinefunction --------------------------------------------------------

Zcombinefunction <- function(X,fac.level,ord.fac,Gorder,
                             indTwo=NULL, indThree=NULL){
    
    
    Zonefunction <-function(Xone,nlevel,ord=FALSE){
        if(ord==FALSE){
            qp <- nlevel*(nlevel-1)/2
        }else if(ord==TRUE){
            qp <- (nlevel-1)
        }
        zero.mat <- matrix(0,nrow=nrow(Xone),ncol=(2*qp))
        Zone <- cbind(Xone, -Xone, zero.mat)
        return(Zone)
    }
    
    
    Ztwofunction <- function(Xtwo,alevel,blevel,ordAB=c(FALSE,FALSE)){
        if(ordAB[1]==FALSE){
            r1 <- (blevel)*(alevel)*(alevel-1)/2
        }else if(ordAB[1]==TRUE){
            r1 <- (blevel)*(alevel-1)
        }
        if(ordAB[2]==FALSE){
            r2 <- (alevel)*(blevel)*(blevel-1)/2
        }else if(ordAB[2]==TRUE){
            r2 <- (alevel)*(blevel-1)
        }    
        zero.mat1 <- matrix(0,nrow=nrow(Xtwo),ncol=(2*r1))
        zero.mat2 <- matrix(0,nrow=nrow(Xtwo),ncol=(2*r2))
        Ztwo <- cbind(Xtwo, -Xtwo, zero.mat1, zero.mat2)
        return(Ztwo)
    }
    
    Zthreefunction <- function(Xthree,alevel,blevel,clevel,ordABC=c(FALSE,FALSE,FALSE)){
        if(ordABC[1]==FALSE){
            r1 <- (blevel*clevel)*(alevel)*(alevel-1)/2
        }else if(ordABC[1]==TRUE){
            r1 <- (blevel*clevel)*(alevel-1)
        }
        if(ordABC[2]==FALSE){
            r2 <- (alevel*clevel)*(blevel)*(blevel-1)/2
        }else if(ordABC[2]==TRUE){
            r2 <- (alevel*clevel)*(blevel-1)
        }
        if(ordABC[3]==FALSE){
            r3 <- (alevel*blevel)*(clevel)*(clevel-1)/2
        }else if(ordABC[3]==TRUE){
            r3 <- (alevel*blevel)*(clevel-1)
        }
        zero.mat1 <- matrix(0,nrow=nrow(Xthree),ncol=(2*r1))
        zero.mat2 <- matrix(0,nrow=nrow(Xthree),ncol=(2*r2))
        zero.mat3 <- matrix(0,nrow=nrow(Xthree),ncol=(2*r3))
        Zthree <- cbind(Xthree, -Xthree, zero.mat1, zero.mat2, zero.mat3)
        return(Zthree)
    }
    ## Input should not have a intercept
    ## Add intercept later
    ## This does not change with ordered categorical variables.
    
    ## fac.level <- c(4,7,4)
    ## ord <- c(TRUE,TRUE,TRUE)
    
    levelIndex <- CreatelevelIndex(fac.level=fac.level,ord.fac=ord.fac,Gorder=Gorder,
                                   indTwo=indTwo, indThree=indThree)
    use.ind <- (levelIndex$plus==1) * (levelIndex$dif==0)
    Index.use <- levelIndex[use.ind==1,]
    Index.use$start <- c(1,cumsum(Index.use$length)[-nrow(Index.use)]+1)
    Index.use$end <-  cumsum(Index.use$length)
    
    Fac.index <- Index.use[,regexpr("Fac",colnames(Index.use))>0]
    
    Z.list <- list()
    
    for(i in 1:nrow(Index.use)){
        if(Index.use$order[i]==1){
            Xone <- X[,Index.use$start[i]:Index.use$end[i]]
            level.main <- fac.level[which(Fac.index[i,]==1)]
            ord.main <- ord.fac[which(Fac.index[i,]==1)]
            Z.list[[i]] <- Zonefunction(Xone=Xone,nlevel=level.main,ord=ord.main)
        }else if(Index.use$order[i]==2){
            Xtwo <- X[,Index.use$start[i]:Index.use$end[i]]
            level.main <- fac.level[which(Fac.index[i,]==1)]
            ord.main <- ord.fac[which(Fac.index[i,]==1)]
            Z.list[[i]] <- Ztwofunction(Xtwo=Xtwo,alevel=level.main[1],
                                        blevel=level.main[2],ordAB=ord.main)
        }else if(Index.use$order[i]==3){
            Xthree <- X[,Index.use$start[i]:Index.use$end[i]]
            level.main <- fac.level[which(Fac.index[i,]==1)]
            ord.main <- ord.fac[which(Fac.index[i,]==1)]
            Z.list[[i]] <- Zthreefunction(Xthree=Xthree,
                                          alevel=level.main[1],blevel=level.main[2],
                                          clevel=level.main[3],
                                          ordABC=ord.main)
        }
    }
    
    Zcombine.temp <- do.call(cbind,Z.list)
    Zcombine <- cbind(1, Zcombine.temp)
    
    return(Zcombine)
}

# Glsei -------------------------------------------------------------------

Glsei <- function (A = NULL, B = NULL, E = NULL, F = NULL, G = NULL, H = NULL, 
                   Wx = NULL, Wa = NULL, type = 2, tol = sqrt(.Machine$double.eps), 
                   tolrank = NULL, fulloutput = FALSE, verbose = TRUE) 
{
    if (is.vector(E) & length(F) == 1) 
        E <- t(E)
    else if (!is.matrix(E) & !is.null(E)) 
        E <- as.matrix(E)
    if (is.vector(A) & length(B) == 1) 
        A <- t(A)
    else if (!is.matrix(A) & !is.null(A)) 
        A <- as.matrix(A)
    if (is.vector(G) & length(H) == 1) 
        G <- t(G)
    else if (!is.matrix(G) & !is.null(G)) 
        G <- as.matrix(G)
    if (!is.matrix(F) & !is.null(F)) 
        F <- as.matrix(F)
    if (!is.matrix(B) & !is.null(B)) 
        B <- as.matrix(B)
    if (!is.matrix(H) & !is.null(H)) 
        H <- as.matrix(H)
    if (is.null(A) && is.null(E)) {
        if (is.null(G)) 
            stop("cannot solve least squares problem - A, E AND G are NULL")
        A <- matrix(data = 0, nrow = 1, ncol = ncol(G))
        B <- 0
    }
    else if (is.null(A)) {
        A <- matrix(data = E[1, ], nrow = 1)
        B <- F[1]
    }
    Neq <- nrow(E)
    Napp <- nrow(A)
    Nx <- ncol(A)
    Nin <- nrow(G)
    if (is.null(Nx)) 
        Nx <- ncol(E)
    if (is.null(Nx)) 
        Nx <- ncol(G)
    if (is.null(Nin)) 
        Nin <- 1
    if (is.null(Neq)) {
        Neq <- 0
        if (verbose & type == 1) 
            warning("No equalities - setting type = 2")
        type = 2
        F <- NULL
    }
    else {
        if (ncol(E) != Nx) 
            stop("cannot solve least squares problem - A and E not compatible")
        if (length(F) != Neq) 
            stop("cannot solve least squares problem - E and F not compatible")
    }
    if (is.null(G)) 
        G <- matrix(data = 0, nrow = 1, ncol = Nx)
    if (is.null(H)) 
        H <- 0
    if (ncol(G) != Nx) 
        stop("cannot solve least squares problem - A and G not compatible")
    if (length(B) != Napp) 
        stop("cannot solve least squares problem - A and B not compatible")
    if (length(H) != Nin) 
        stop("cannot solve least squares problem - G and H not compatible")
    if (!is.null(Wa)) {
        if (length(Wa) != Napp) 
            stop("cannot solve least squares problem - Wa should have length = number of rows of A")
        A <- A * Wa
        B <- B * Wa
    }
    Tol <- tol
    if (is.null(Tol)) 
        Tol <- sqrt(.Machine$double.eps)
    IsError <- FALSE
    if (type == 2) {
        if (!is.null(Wx)) 
            stop("cannot solve least squares problem - weights not implemented for type 2")
        if (!is.null(Wa)) 
            stop("cannot solve least squares problem - weights not implemented for type 2")
        ## library("quadprog")
        dvec <- crossprod(A, B)
        Dmat <- crossprod(A, A)
        Pass <- 0
        tol.d <- tol
        while(Pass==0){
            diag(Dmat) <- diag(Dmat) + tol.d
            Amat <- t(rbind(E, G))
            bvec <- c(F, H)
            sol <- solve.QP(Dmat, dvec, Amat, bvec, meq = Neq)
            sol$IsError <- FALSE
            sol$X <- sol$solution
            if(any(is.na(sol$X))==TRUE){
                Pass <- 0
                tol.d <- 5*tol.d
            }else{
                Pass <- 1
            }
        }
    }
    else stop("cannot solve least squares problem - type unknown")
    X <- sol$X
    X[which(abs(X) < Tol)] <- 0
    if (any(is.infinite(X))) {
        residual <- Inf
        solution <- Inf
    }
    else {
        residual <- 0
        if (Nin > 0) {
            ineq <- G %*% X - H
            residual <- residual - sum(ineq[ineq < 0])
        }
        if (Neq > 0) 
            residual <- residual + sum(abs(E %*% X - F))
        if (residual > Tol) 
            sol$IsError <- TRUE
        solution <- 0
        if (Napp > 0) 
            solution <- sum((A %*% X - B)^2)
    }
    xnames <- colnames(A)
    if (is.null(xnames)) 
        xnames <- colnames(E)
    if (is.null(xnames)) 
        xnames <- colnames(G)
    names(X) <- xnames
    res <- list(X = X, residualNorm = residual, solutionNorm = solution, 
                IsError = sol$IsError, type = "lsei")
    return(res)
}


# Collapsing --------------------------------------------------------------

Collapsing <- function(fit){
    ## Make the difference in three-way interactions
    D3function <- function(alevel,blevel,clevel,type="A",ordABC=c(FALSE,FALSE,FALSE)){
        if(type=="A"){
            D3 <- D3Afunction(alevel=alevel,blevel=blevel,clevel=clevel,ord=ordABC[1])
        }else if(type=="B"){
            D3 <- D3Bfunction(alevel=alevel,blevel=blevel,clevel=clevel,ord=ordABC[2])
        }else if(type=="C"){
            D3 <- D3Cfunction(alevel=alevel,blevel=blevel,clevel=clevel,ord=ordABC[3])
        }
        return(D3)
    }
    
    D3Afunction <- function(alevel,blevel,clevel,ord=FALSE){
        D3a <- matrix(0,nrow=0,ncol=(alevel*blevel*clevel))  
        if(ord==FALSE){
            for(i in 1:(alevel-1)){
                for(j  in (i+1):alevel){
                    base.v <- rep(0,alevel)
                    base.v[i] <- -1
                    base.v[j] <- 1
                    base.mat <- diag((blevel*clevel)) %x% t(base.v)
                    D3a <- rbind(D3a,base.mat)
                }
            }
        }else if(ord==TRUE){
            for(i in 1:(alevel-1)){
                base.v <- rep(0,alevel)
                base.v[i] <- -1
                base.v[(i+1)] <- 1
                base.mat <- diag((blevel*clevel)) %x% t(base.v)
                D3a <- rbind(D3a,base.mat)
            }        
        }
        return(D3a)
    }
    
    ## (B): Make the difference in two-way interactions 
    D3Bfunction <- function(alevel,blevel,clevel,ord=FALSE){    
        if(ord==FALSE){
            D3b <- matrix(0,nrow=0,ncol=(alevel*blevel*clevel))
            for(i in 1:(blevel-1)){     
                left.temp <- matrix(0,ncol=((i-1)*(alevel)),nrow=(alevel))
                left.main <- -diag((alevel))
                left.final <- cbind(left.temp,left.main)
                for(j  in (i+1):blevel){
                    middle.temp <-  matrix(0,ncol=((j-i-1)*(alevel)),nrow=(alevel))
                    right.main <- diag((alevel))
                    right.temp <- matrix(0,ncol=((blevel-j)*(alevel)),nrow=(alevel))
                    right.final <- cbind(middle.temp,right.main,right.temp)
                    D2b <- cbind(left.final,right.final) ## The base
                    D3base <- diag(clevel) %x% D2b
                    D3b <- rbind(D3b, D3base)
                }
            }
        }else if (ord==TRUE){
            D3b <- matrix(0,nrow=0,ncol=(alevel*blevel*clevel))
            for(i in 1:(blevel-1)){     
                left <- matrix(0,ncol=((i-1)*(alevel)),nrow=(alevel))
                Main <- cbind(-diag((alevel)),diag((alevel)))
                right <- matrix(0,ncol=((blevel-(i+1))*(alevel)),nrow=(alevel))
                D2b <- cbind(left,Main,right)
                D3base <- diag(clevel) %x% D2b
                D3b <- rbind(D3b, D3base)
            }
        }
        return(D3b)
        ## psy_B(2,3), psy_B(2,4), psy_B(2,5) 
    }
    
    D3Cfunction <- function(alevel,blevel,clevel,ord=FALSE){
        if(ord==FALSE){
            D3c <- matrix(0,nrow=0,ncol=(alevel*blevel*clevel))
            for(i in 1:(clevel-1)){     
                left.temp <- matrix(0,ncol=((i-1)*(alevel*blevel)),nrow=(alevel*blevel))
                left.main <- -diag((alevel*blevel))
                left.final <- cbind(left.temp,left.main)
                for(j  in (i+1):clevel){
                    middle.temp <-  matrix(0,ncol=((j-i-1)*(alevel*blevel)),nrow=(alevel*blevel))
                    right.main <- diag((alevel*blevel))
                    right.temp <- matrix(0,ncol=((clevel-j)*(alevel*blevel)),
                                         nrow=(alevel*blevel))
                    right.final <- cbind(middle.temp,right.main,right.temp)
                    D3c <- rbind(D3c,cbind(left.final,right.final))
                }
            }
        }else if (ord==TRUE){
            Left <- cbind(-diag((alevel*blevel)*(clevel-1)),
                          matrix(0,ncol=(alevel*blevel),nrow=(alevel*blevel)*(clevel-1)))
            Right <- cbind(matrix(0,ncol=(alevel*blevel),nrow=(alevel*blevel)*(clevel-1)),
                           diag((alevel*blevel)*(clevel-1)))
            D3c <- Left + Right
        }
        return(D3c)
    }
    
    
    ## From this version, I use D2 rather than D3 
    ## Make the difference in two-way interactions 
    D2function <- function(alevel,blevel,type,ordAB){
        if(type=="A"){
            D2 <- D2Afunction(alevel=alevel,blevel=blevel,ord=ordAB[1])
        }else if(type=="B"){        
            D2 <- D2Bfunction(alevel=alevel,blevel=blevel,ord=ordAB[2])
        }
        return(D2)
        ## D2 matrix
        ## nrow = differences between two-way interactions fixing the sublevel
        ## ncol = alevel*blevel (coefficients for two-ways)
    }
    
    ## (A): Make the difference in two-way interactions 
    D2Afunction <- function(alevel,blevel,ord=FALSE){
        D2a <- matrix(0,nrow=0,ncol=(alevel*blevel))
        if(ord==FALSE){
            for(i in 1:(alevel-1)){
                for(j  in (i+1):alevel){
                    base.v <- rep(0,alevel)
                    base.v[i] <- -1
                    base.v[j] <- 1
                    base.mat <- diag(blevel) %x% t(base.v)
                    D2a <- rbind(D2a,base.mat)
                }
            }
        }else if(ord==TRUE){
            for(i in 1:(alevel-1)){
                base.v <- rep(0,alevel)
                base.v[i] <- -1
                base.v[(i+1)] <- 1
                base.mat <- diag(blevel) %x% t(base.v)
                D2a <- rbind(D2a,base.mat)
            }        
        }
        return(D2a)
    }
    
    ## (B): Make the difference in two-way interactions 
    D2Bfunction <- function(alevel,blevel,ord=FALSE){    
        if(ord==FALSE){
            D2b <- matrix(0,nrow=0,ncol=(alevel*blevel))
            for(i in 1:(blevel-1)){     
                left.temp <- matrix(0,ncol=((i-1)*(alevel)),nrow=(alevel))
                left.main <- -diag((alevel))
                left.final <- cbind(left.temp,left.main)
                for(j  in (i+1):blevel){
                    middle.temp <-  matrix(0,ncol=((j-i-1)*(alevel)),nrow=(alevel))
                    right.main <- diag((alevel))
                    right.temp <- matrix(0,ncol=((blevel-j)*(alevel)),nrow=(alevel))
                    right.final <- cbind(middle.temp,right.main,right.temp)
                    D2b <- rbind(D2b,cbind(left.final,right.final))
                }
            }
        }else if (ord==TRUE){
            Left <- cbind(-diag((alevel)*(blevel-1)),matrix(0,ncol=alevel,nrow=(alevel)*(blevel-1)))
            Right <- cbind(matrix(0,ncol=alevel,nrow=(alevel)*(blevel-1)), diag((alevel)*(blevel-1)))
            D2b <- Left + Right
        }
        return(D2b)
        ## psy_B(2,3), psy_B(2,4), psy_B(2,5) 
    }
    
    ## (Main): Making the difference between all levels
    D1function <- function(nlevel,ord=FALSE){
        if(ord==FALSE){
            qp = (nlevel)*(nlevel-1)/2
            D1.mat <- matrix(0,ncol=(nlevel),nrow=0)    
            if (qp==1) D1.mat <- matrix(c(-1,1),ncol=2,nrow=1)
            if (qp>1){
                for(i in 1 : (nlevel-1)){
                    w.diag <- diag((nlevel-i))
                    left.v <- c(rep(0,(i-1)),-1)
                    left.mat <- matrix(rep(left.v,(nlevel-i)),nrow=(nlevel-i),byrow=TRUE)
                    w.mat <- cbind(left.mat,w.diag)
                    D1.mat <- rbind(D1.mat,w.mat)
                }
            }
        }else if(ord==TRUE){
            if(nlevel==2){
                D1.mat <- matrix(c(-1,1),ncol=2,nrow=1)
            }else{
                D1.mat <- cbind(0,diag(nlevel-1)) + cbind(-diag(nlevel-1),0)
            }
            
        }    
        return(D1.mat)
        ## nrow = differences between main effects
        ## ncol = mainlevel (coefficients for the main factor)
    }
    
    fac.level <- fit$fac.level
    ord.fac  <- fit$ord.fac
    AME <- fit$AME
    AMIE2 <- fit$AMIE2
    AMIE3 <- fit$AMIE3
    eps <- fit$eps   
    fac.name <- all.vars(fit$formula)[-1]
    n.fac <- length(fac.level)
    indTwo <- fit$indTwo
    indThree <- fit$indThree
    Gorder <- fit$Gorder
    
    levelIndex <- CreatelevelIndex(fac.level=fac.level,ord.fac=ord.fac,Gorder=Gorder,
                                   indTwo=indTwo, indThree=indThree)
    use.ind <-  (levelIndex$plus==1)*(levelIndex$dif==0)
    Index.use <- levelIndex[use.ind==1,]
    Fac.index <- levelIndex[,regexpr("Fac",colnames(levelIndex))>0]
    Fac.Ind.use <- Fac.index[use.ind==1,]
    
    order.f <- attr(terms(fit$formula, data=fit$data), "order")
    Fac.Ind.use1 <- Fac.Ind.use[order.f==1,]
    if(any(order.f==2)) Fac.Ind.use2 <- Fac.Ind.use[order.f==2,]
    if(any(order.f==3)) Fac.Ind.use3 <- Fac.Ind.use[order.f==3,]
    
    
    ## Do Collapsing for each factor
    Collapse <- list()
    for(z in 1:n.fac){
        ## First Order 
        ## type.ind <- z    
        MainD1 <- D1function(nlevel=fac.level[z],ord=ord.fac[z])
        MainDif <- t(MainD1 %*% AME[[z]])
        
        ## Two-ways
        onlyOne <- sum(Fac.Ind.use[,z])==1
        if(any(order.f==2)==TRUE) yesTwo <-  sum(Fac.Ind.use2[,z]==1)>0 else yesTwo <- FALSE
        if(any(order.f==3)==TRUE) yesThree <- sum(Fac.Ind.use3[,z]==1)>0 else yesThree <- FALSE
        if(onlyOne==TRUE){
            ## No Interaction
            Dif <- abs(MainDif)
            Collapse[[z]] <- apply(Dif <= eps, 2, all)
            ## print(Dif)
        }else if(yesTwo==TRUE){
            
            int2.index <- which(Fac.Ind.use2[,z]==1)
            D2.mat <- matrix(NA, ncol=ncol(MainDif), nrow=0)
            for(w in int2.index){
                ## Setup the coefficients
                AMIE2.u <- AMIE2[[w]]
                
                ## Construct D matrix
                Fac1 <- min(which(Fac.Ind.use2[w,]==1))
                Fac2 <- max(which(Fac.Ind.use2[w,]==1))
                if(z == Fac1) type.d2 <- "A"
                if(z == Fac2) type.d2 <- "B"
                D2  <- D2function(alevel=fac.level[Fac1],blevel=fac.level[Fac2],
                                  type=type.d2, ordAB=ord.fac[c(Fac1,Fac2)])
                D2.mat0 <- D2 %*% AMIE2.u
                D2.mat01  <- matrix(D2.mat0, ncol=ncol(MainDif))
                D2.mat <- rbind(D2.mat, D2.mat01)
            }
            D2.mat <- abs(D2.mat)
            if(yesThree==FALSE){                
                Dif <- abs(rbind(MainDif, D2.mat))
                Collapse[[z]] <- apply(Dif <= eps, 2, all)
            }else if(yesThree==TRUE){
                int3.index <- which(Fac.Ind.use3[,z]==1)
                D3.mat <- matrix(NA, ncol=ncol(MainDif), nrow=0)
                
                for(w in int3.index){
                    ## Setup the coefficients
                    AMIE3.u <- AMIE3[[w]]
                    
                    ## Construct D matrix
                    Fac1 <- which(Fac.Ind.use3[w,]==1)[1]
                    Fac2 <- which(Fac.Ind.use3[w,]==1)[2]
                    Fac3 <- which(Fac.Ind.use3[w,]==1)[3]
                    if(z == Fac1) type.d3 <- "A"
                    if(z == Fac2) type.d3 <- "B"
                    if(z == Fac3) type.d3 <- "C"
                    D3  <- D3function(alevel=fac.level[Fac1], blevel=fac.level[Fac2], clevel=fac.level[Fac3],
                                      type=type.d3, ordAB=ord.fac[c(Fac1,Fac2, Fac3)])
                    D3.mat0 <- D3 %*% c(AMIE3.u)
                    D3.mat01  <- matrix(D3.mat0, ncol=ncol(MainDif))
                    D3.mat <- rbind(D3.mat, D3.mat01)
                }
                D3.mat <- abs(D3.mat)         
                Dif <- abs(rbind(MainDif, D2.mat, D3.mat))
                Collapse[[z]] <- apply(Dif <= eps, 2, all)
            }
        }
    }
    
    ## I got Collapsing Index, then decide which levels will be collapsed.
    collapse.level <- list()
    for(z in 1:n.fac){
        adj <- matrix(0, ncol=fac.level[z], nrow=fac.level[z])
        if(ord.fac[z]==TRUE){
            for(i in 1:length(Collapse[[z]])){
                adj[i, (i+1)] <- as.numeric(Collapse[[z]][i])
            }
            adj <- adj + t(adj)
        }else if(ord.fac[z]==FALSE){
            ref <- combn(seq(1:fac.level[z]),2)
            for(i in 1:length(Collapse[[z]])){
                adj[ref[1,i], ref[2,i]] <- as.numeric(Collapse[[z]][i])
            }
            adj <- adj + t(adj)
        }        
        g <- graph_from_adjacency_matrix(adj, mode="undirected")
        collapse.level[[z]] <- components(g)$membership
    }
    names(collapse.level) <- fac.name
    
    ## Combine two weights. 
    output <- list("Collapse.Index"=Collapse, "collapse.level"=collapse.level)
    return(output)
}
